4DNucleome · Czaki · Mar 13, 2025 · Mar 12, 2025 · Mar 12, 2025 · Mar 13, 2025
diff --git a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp
@@ -112,19 +112,27 @@ inline int64_t double_as_hex(double d) {
   return result;
 }
 
+enum Orientation {
+  COLLINEAR = 0,
+  CLOCKWISE = 1,
+  COUNTERCLOCKWISE = 2,
+};
+
 /**
  * Determines the orientation of the triplet (p, q, r).
  *
  * @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.
+ * @return Orientation::COLLINEAR if p, q and r are collinear.
+ *         Orientation::CLOCKWISE if the triplet (p, q, r) is in a clockwise
+ *         orientation.
+ *         Orientation::COUNTERCLOCKWISE if the triplet (p, q, r) is in a
+ *         counterclockwise orientation.
  */
-inline int _orientation(const point::Point &p, const point::Point &q,
-                        const point::Point &r) {
+inline Orientation _orientation(const point::Point &p, const point::Point &q,
+                                const point::Point &r) {
   double val1 = ((q.y - p.y) * (r.x - q.x));
   double val2 = ((r.y - q.y) * (q.x - p.x));
   // This commented code if for debugging purposes of differences between macOS
@@ -142,8 +150,8 @@ inline int _orientation(const point::Point &p, const point::Point &q,
   // }
   // Instead of using classical equation, we need to use two variables
   // to handle problem with strange behavior on macOS.
-  if (val1 == val2) return 0;
-  return (val1 > val2) ? 1 : 2;
+  if (val1 == val2) return Orientation::COLLINEAR;
+  return (val1 > val2) ? Orientation::CLOCKWISE : Orientation::COUNTERCLOCKWISE;
 }
 
 /**
@@ -163,17 +171,17 @@ inline bool _do_intersect(const point::Segment &s1, const point::Segment &s2) {
   const point::Point &p2 = s2.bottom;
   const point::Point &q2 = s2.top;
 
-  int o1 = _orientation(p1, q1, p2);
-  int o2 = _orientation(p1, q1, q2);
-  int o3 = _orientation(p2, q2, p1);
-  int o4 = _orientation(p2, q2, q1);
+  Orientation o1 = _orientation(p1, q1, p2);
+  Orientation o2 = _orientation(p1, q1, q2);
+  Orientation o3 = _orientation(p2, q2, p1);
+  Orientation o4 = _orientation(p2, q2, q1);
 
   if (o1 != o2 && o3 != o4) return true;
 
-  if (o1 == 0 && _on_segment(p1, p2, q1)) return true;
-  if (o2 == 0 && _on_segment(p1, q2, q1)) return true;
-  if (o3 == 0 && _on_segment(p2, p1, q2)) return true;
-  if (o4 == 0 && _on_segment(p2, q1, q2)) return true;
+  if (o1 == Orientation::COLLINEAR && _on_segment(p1, p2, q1)) return true;
+  if (o2 == Orientation::COLLINEAR && _on_segment(p1, q2, q1)) return true;
+  if (o3 == Orientation::COLLINEAR && _on_segment(p2, p1, q2)) return true;
+  if (o4 == Orientation::COLLINEAR && _on_segment(p2, q1, q2)) return true;
 
   return false;
 }

diff --git a/src/PartSegCore_compiled_backend/triangulation/point.hpp b/src/PartSegCore_compiled_backend/triangulation/point.hpp
@@ -208,6 +208,19 @@ struct Segment {
     }
   };
 };
+Point centroid(const std::vector<Point> &point_list) {
+  if (point_list.empty()) {
+    return {0, 0};
+  }
+  Point res(0, 0);
+  for (auto &point : point_list) {
+    res.x += point.x;
+    res.y += point.y;
+  }
+  res.x /= float(point_list.size());
+  res.y /= float(point_list.size());
+  return res;
+}
 }  // namespace partsegcore::point
 
 // overload of hash function for

diff --git a/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp
@@ -2,6 +2,7 @@
 #define PARTSEGCORE_TRIANGULATE_H
 
 #include <algorithm>
+#include <cmath>
 #include <map>
 #include <memory>  // memory header is required on linux, and not on macos
 #include <set>
@@ -602,6 +603,36 @@ struct MonotonePolygonBuilder {
   };
 };
 
+/*
+ * Check if the polygon, that all angles have the same orientation
+ * do not have self-intersections
+ *
+ * @param begin Iterator to the first point of the polygon
+ * @param end Iterator to the end of the polygon
+ * @param centroid Centroid of the polygon
+ *
+ * @return True if the polygon is simple, false otherwise
+ */
+template <typename Iterator>
+bool is_simple_polygon(Iterator begin, Iterator end, point::Point centroid) {
+  double start_angle = std::atan2(begin->y - centroid.y, begin->x - centroid.x);
+  double prev_angle = 0;
+  begin++;
+  for (; begin != end; begin++) {
+    double angle =
+        std::atan2(begin->y - centroid.y, begin->x - centroid.x) - start_angle;
+    if (angle < 0) {
+      angle += 2 * M_PI;
+    }
+    if (angle < prev_angle) {
+      return false;
+    } else {
+      prev_angle = angle;
+    }
+  }
+  return true;
+}
+
 /**
  * Checks if a given polygon is convex.
  *
@@ -615,26 +646,47 @@ struct MonotonePolygonBuilder {
  * @return True if the polygon is convex, false otherwise.
  */
 inline bool _is_convex(const std::vector<point::Point> &polygon) {
-  int orientation = 0;
-  int triangle_orientation;
-  for (std::size_t i = 0; i < polygon.size() - 2; i++) {
-    triangle_orientation =
-        intersection::_orientation(polygon[i], polygon[i + 1], polygon[i + 2]);
-    if (triangle_orientation == 0) continue;
-    if (orientation == 0)
+  if (polygon.size() < 3) return false;
+  if (polygon.size() == 3) return true;
+  intersection::Orientation orientation = intersection::Orientation::COLLINEAR;
+  intersection::Orientation triangle_orientation;
+  size_t idx = 0;
+  for (; idx < polygon.size() - 2; idx++) {
+    triangle_orientation = intersection::_orientation(
+        polygon[idx], polygon[(idx + 1)], polygon[(idx + 2)]);
+    if (triangle_orientation != intersection::Orientation::COLLINEAR) {
       orientation = triangle_orientation;
-    else if (orientation != triangle_orientation)
+      break;
+    }
+  }
+  if (orientation == intersection::Orientation::COLLINEAR) {
+    return false;
+  }
+  for (; idx < polygon.size() - 2; idx++) {
+    triangle_orientation = intersection::_orientation(
+        polygon[idx], polygon[(idx + 1)], polygon[(idx + 2)]);
+    if (triangle_orientation != 0 && triangle_orientation != orientation) {
       return false;
+    }
   }
   triangle_orientation = intersection::_orientation(
       polygon[polygon.size() - 2], polygon[polygon.size() - 1], polygon[0]);
-  if (triangle_orientation != 0 && triangle_orientation != orientation)
+  if (triangle_orientation != 0 && triangle_orientation != orientation) {
     return false;
+  }
   triangle_orientation = intersection::_orientation(polygon[polygon.size() - 1],
                                                     polygon[0], polygon[1]);
-  if (triangle_orientation != 0 && triangle_orientation != orientation)
+  if (triangle_orientation != 0 && triangle_orientation != orientation) {
     return false;
-  return true;
+  }
+
+  point::Point centroid = point::centroid(polygon);
+
+  if (orientation == intersection::Orientation::COUNTERCLOCKWISE) {
+    return is_simple_polygon(polygon.begin(), polygon.end(), centroid);
+  } else {
+    return is_simple_polygon(polygon.rbegin(), polygon.rend(), centroid);
+  }
 }
 
 /**

diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py
@@ -208,6 +208,7 @@ def test_triangle_convex_polygon():
     ],
 )
 def test_triangulate_polygon_py_convex(polygon, expected):
+    assert is_convex(polygon)
     assert triangulate_polygon_py(polygon)[0] == expected
 
 
@@ -691,3 +692,55 @@ def test_splitting_edges(split_edges, triangles):
     )
     triangles_ = triangulate_polygon_with_edge_numpy_li([polygon], split_edges=split_edges)[1][2]
     assert len(triangles_) == triangles
+
+
+def generate_regular_polygon(n: int, reverse: bool, radius: int = 1) -> np.ndarray:
+    angles = np.linspace(0, 2 * np.pi, n, endpoint=False)
+    if reverse:
+        angles = angles[::-1]
+    return np.column_stack((radius * np.cos(angles), radius * np.sin(angles)))
+
+
+def generate_self_intersecting_polygon(n: int, reverse: bool, radius: int = 1) -> np.ndarray:
+    """Generate self-intersecting polygon with n vertices
+
+    The polygon is generated by doubling the angle range of a regular polygon
+    """
+    assert n % 2 == 1, 'an odd number is required to generate a self-intersecting polygon'
+    angles = np.linspace(0, 4 * np.pi, n, endpoint=False)
+    if reverse:
+        angles = angles[::-1]
+    return np.column_stack((radius * np.cos(angles), radius * np.sin(angles)))
+
+
+def rotation_matrix(angle: float) -> np.ndarray:
+    """Create a 2D rotation matrix for the given angle in degrees."""
+    return np.array(
+        [
+            [np.cos(np.radians(angle)), -np.sin(np.radians(angle))],
+            [np.sin(np.radians(angle)), np.cos(np.radians(angle))],
+        ]
+    )
+
+
+ANGLES = [0, 5, 75, 95, 355]
+
+
+@pytest.mark.parametrize('angle', ANGLES, ids=str)
+@pytest.mark.parametrize('n_vertex', [5, 7, 19])
+@pytest.mark.parametrize('reverse', [False, True])
+def test_is_convex_self_intersection(angle, n_vertex, reverse):
+    p = generate_self_intersecting_polygon(n_vertex, reverse)
+    rot = rotation_matrix(angle)
+    data = np.dot(p, rot)
+    assert not is_convex(data)
+
+
+@pytest.mark.parametrize('angle', ANGLES, ids=str)
+@pytest.mark.parametrize('n_vertex', [3, 4, 7, 12, 15, 20])
+@pytest.mark.parametrize('reverse', [False, True])
+def test_is_convex_regular_polygon(angle, n_vertex, reverse):
+    poly = generate_regular_polygon(n_vertex, reverse=reverse)
+    rot = rotation_matrix(angle)
+    rotated_poly = np.dot(poly, rot)
+    assert is_convex(rotated_poly)