add Cython implementation of is_point_in_polygon_2d()

Author: mozman

profiling/is_point_in_polygon.py

@@ -0,0 +1,69 @@
+#  Copyright (c) 2024, Manfred Moitzi
+#  License: MIT License
+import time
+
+from ezdxf.math import Vec2
+from ezdxf.math._construct import is_point_in_polygon_2d
+from ezdxf.acc.construct import is_point_in_polygon_2d as is_point_in_polygon_cy
+
+VERTICES = [(1, 1), (5, 1), (5, 3), (3, 3), (3, 5), (5, 5), (5, 7), (1, 7)]
+SHAPE_PY = Vec2.list(VERTICES)
+
+POINTS_INSIDE = Vec2.list(
+    [(2, 2), (3, 2), (4, 2), (2, 3), (2, 4), (2, 5), (2, 6), (3, 6), (4, 6)]
+)
+
+POINTS_OUTSIDE = Vec2.list(
+    [
+        (0, 0),  # 0
+        (2, 0),  # 1
+        (4, 0),  # 2
+        (6, 0),  # 3
+        (0, 2),  # 4
+        (6, 2),  # 5
+        (0, 4),  # 6
+        (4, 4),  # 7
+        (6, 2),  # 8
+        (0, 6),  # 9
+        (6, 6),  # a
+        (0, 8),  # b
+        (2, 8),  # c
+        (4, 8),  # d
+        (6, 8),  # e
+    ]
+)
+
+
+ROUNDS = 2000
+
+
+def python_version():
+    for point in POINTS_INSIDE + POINTS_OUTSIDE:
+        is_point_in_polygon_2d(point, SHAPE_PY)
+
+
+def cython_version():
+    for point in POINTS_INSIDE + POINTS_OUTSIDE:
+        is_point_in_polygon_cy(point, SHAPE_PY)
+
+
+def profile(func) -> float:
+    t0 = time.perf_counter()
+    for _ in range(ROUNDS):
+        func()
+    t1 = time.perf_counter()
+    return t1 - t0
+
+
+def main():
+    py_t = profile(python_version)
+    print(f"python version: {py_t:.3f}s")
+
+    # Numpy/Cython version was just 10x faster, this Cython version is 23x faster!
+    cy_t = profile(cython_version)
+    print(f"cython version: {cy_t:.3f}s")
+    print(f"ratio python/cython: {py_t/cy_t:.3f}")
+
+
+if __name__ == "__main__":
+    main()

src/ezdxf/acc/construct.pyx

@@ -1,6 +1,6 @@
 # cython: language_level=3
 # distutils: language = c++
-# Copyright (c) 2020-2023, Manfred Moitzi
+# Copyright (c) 2020-2024, Manfred Moitzi
 # License: MIT License
 # type: ignore -- pylance sucks at type-checking cython files
 from typing import Iterable, TYPE_CHECKING, Sequence, Optional, Tuple
@@ -199,3 +199,80 @@ def arc_angle_span_rad(double start, double end) -> float:
     if end < start:
         end += M_TAU
     return end - start
+
+
+def is_point_in_polygon_2d(
+    point: Vec2, polygon: list[Vec2], double abs_tol=TOLERANCE
+) -> int:
+    """
+    Test if `point` is inside `polygon`.  Returns +1 for inside, 0 for on the 
+    boundary and  -1 for outside.
+
+    Supports convex and concave polygons with clockwise or counter-clockwise oriented
+    polygon vertices.  Does not raise an exception for degenerated polygons.
+
+
+    Args:
+        point: 2D point to test as :class:`Vec2`
+        polygon: list of 2D points as :class:`Vec2`
+        abs_tol: tolerance for distance check
+
+    Returns:
+        +1 for inside, 0 for on the boundary, -1 for outside
+
+    """
+    # Source: http://www.faqs.org/faqs/graphics/algorithms-faq/
+    # Subject 2.03: How do I find if a point lies within a polygon?
+    # Numpy version was just 10x faster, this version is 23x faster than the Python 
+    # version!
+    cdef double  a, b, c, d, x, y, x1, y1, x2, y2
+    cdef list vertices = polygon
+    cdef Vec2 p1, p2
+    cdef int size, last, i
+    cdef bint inside = 0
+
+    size = len(vertices)
+    if size < 3:  # empty polygon
+        return -1
+    last = size - 1
+    p1 = <Vec2> vertices[0]
+    p2 = <Vec2> vertices[last]
+
+    if v2_isclose(p1, p2, REL_TOL, ABS_TOL):  # open polygon
+        size -= 1
+        last -= 1
+    if size < 3:
+        return -1
+
+    x = point.x
+    y = point.y
+    p1 = <Vec2> vertices[last]
+    x1 = p1.x
+    y1 = p1.y
+
+    for i in range(size):
+        p2 = <Vec2> vertices[i]
+        x2 = p2.x
+        y2 = p2.y
+
+        # is point on polygon boundary line:
+        # is point in x-range of line
+        a, b = (x2, x1) if x2 < x1 else (x1, x2)
+        if a <= x <= b:
+            # is point in y-range of line
+            c, d = (y2, y1) if y2 < y1 else (y1, y2)
+            if (c <= y <= d) and fabs(
+                (y2 - y1) * x - (x2 - x1) * y + (x2 * y1 - y2 * x1)
+            ) <= abs_tol:
+                return 0  # on boundary line
+        if ((y1 <= y < y2) or (y2 <= y < y1)) and (
+            x < (x2 - x1) * (y - y1) / (y2 - y1) + x1
+        ):
+            inside = not inside
+        x1 = x2
+        y1 = y2
+    if inside:
+        return 1  # inside polygon
+    else:
+        return -1  # outside polygon
+    

src/ezdxf/math/_construct.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2011-2022, Manfred Moitzi
+# Copyright (c) 2011-2024, Manfred Moitzi
 # License: MIT License
 # These are the pure Python implementations of the Cython accelerated
 # construction tools: ezdxf/acc/construct.pyx
@@ -216,3 +216,61 @@ def arc_angle_span_rad(start: float, end: float) -> float:
     if end < start:
         end += tau
     return end - start
+
+
+def is_point_in_polygon_2d(
+    point: Vec2, polygon: list[Vec2], abs_tol=TOLERANCE
+) -> int:
+    """
+    Test if `point` is inside `polygon`.  Returns +1 for inside, 0 for on the 
+    boundary and  -1 for outside.
+
+    Supports convex and concave polygons with clockwise or counter-clockwise oriented
+    polygon vertices.  Does not raise an exception for degenerated polygons.
+
+
+    Args:
+        point: 2D point to test as :class:`Vec2`
+        polygon: list of 2D points as :class:`Vec2`
+        abs_tol: tolerance for distance check
+
+    Returns:
+        +1 for inside, 0 for on the boundary, -1 for outside
+
+    """
+    # Source: http://www.faqs.org/faqs/graphics/algorithms-faq/
+    # Subject 2.03: How do I find if a point lies within a polygon?
+    # polygon: the Cython implementation needs a list as input to be fast!
+    assert isinstance(polygon, list)
+    if len(polygon) < 3: # empty polygon
+        return -1
+
+    if polygon[0].isclose(polygon[-1]):  # open polygon is required
+        polygon = polygon[:-1]
+    if len(polygon) < 3:
+        return -1
+    x = point.x
+    y = point.y
+    inside = False
+    x1, y1 = polygon[-1]
+    for x2, y2 in polygon:
+        # is point on polygon boundary line:
+        # is point in x-range of line
+        a, b = (x2, x1) if x2 < x1 else (x1, x2)
+        if a <= x <= b:
+            # is point in y-range of line
+            c, d = (y2, y1) if y2 < y1 else (y1, y2)
+            if (c <= y <= d) and abs(
+                (y2 - y1) * x - (x2 - x1) * y + (x2 * y1 - y2 * x1)
+            ) <= abs_tol:
+                return 0  # on boundary line
+        if ((y1 <= y < y2) or (y2 <= y < y1)) and (
+            x < (x2 - x1) * (y - y1) / (y2 - y1) + x1
+        ):
+            inside = not inside
+        x1 = x2
+        y1 = y2
+    if inside:
+        return 1  # inside polygon
+    else:
+        return -1  # outside polygon

src/ezdxf/math/_ctypes.py

@@ -29,6 +29,7 @@
     "intersection_ray_ray_3d",
     "arc_angle_span_deg",
     "arc_angle_span_rad",
+    "is_point_in_polygon_2d",
 ]
 # Import of Python or Cython implementations:
 if USE_C_EXT:
@@ -57,6 +58,7 @@
         intersection_ray_ray_3d,
         arc_angle_span_deg,
         arc_angle_span_rad,
+        is_point_in_polygon_2d,
     )
 else:
     from ._vector import (
@@ -84,6 +86,7 @@
         intersection_ray_ray_3d,
         arc_angle_span_deg,
         arc_angle_span_rad,
+        is_point_in_polygon_2d,
     )
 
 # Early required type aliases

src/ezdxf/math/construct2d.py

@@ -31,7 +31,6 @@
     "is_convex_polygon_2d",
     "is_axes_aligned_rectangle_2d",
     "is_point_on_line_2d",
-    "is_point_in_polygon_2d",
     "is_point_left_of_line",
     "point_to_line_relation",
     "enclosing_angles",
@@ -241,62 +240,8 @@ def distance_point_line_2d(point: Vec2, start: Vec2, end: Vec2) -> float:
     return math.fabs((start - point).det(end - point)) / (end - start).magnitude
 
 
-# Candidate for a faster Cython implementation:
-# is also used for testing 3D ray and line intersection with polygon
-def is_point_in_polygon_2d(
-    point: Vec2, polygon: Sequence[Vec2], abs_tol=TOLERANCE
-) -> int:
-    """Test if `point` is inside `polygon`. Returns ``-1`` (for outside) if the
-    polygon is degenerated, no exception will be raised.
-
-    Supports convex and concave polygons with clockwise or counter-clockwise oriented
-    polygon vertices.
-
-    Args:
-        point: 2D point to test as :class:`Vec2`
-        polygon: sequence of 2D points as :class:`Vec2`
-        abs_tol: tolerance for distance check
-
-    Returns:
-        ``+1`` for inside, ``0`` for on boundary line, ``-1`` for outside
-
-    """
-    # Source: http://www.faqs.org/faqs/graphics/algorithms-faq/
-    # Subject 2.03: How do I find if a point lies within a polygon?
-    # TODO: Cython implementation
-    if not polygon:  # empty polygon
-        return -1
-
-    if polygon[0].isclose(polygon[-1]):  # open polygon is required
-        polygon = polygon[:-1]
-    if len(polygon) < 3:
-        return -1
-    x = point.x
-    y = point.y
-    inside = False
-    x1, y1 = polygon[-1]
-    for x2, y2 in polygon:
-        # is point on polygon boundary line:
-        # is point in x-range of line
-        a, b = (x2, x1) if x2 < x1 else (x1, x2)
-        if a <= x <= b:
-            # is point in y-range of line
-            c, d = (y2, y1) if y2 < y1 else (y1, y2)
-            if (c <= y <= d) and math.fabs(
-                (y2 - y1) * x - (x2 - x1) * y + (x2 * y1 - y2 * x1)
-            ) <= abs_tol:
-                return 0  # on boundary line
-        if ((y1 <= y < y2) or (y2 <= y < y1)) and (
-            x < (x2 - x1) * (y - y1) / (y2 - y1) + x1
-        ):
-            inside = not inside
-        x1 = x2
-        y1 = y2
-    if inside:
-        return 1  # inside polygon
-    else:
-        return -1  # outside polygon
 
+    
 
 def circle_radius_3p(a: Vec3, b: Vec3, c: Vec3) -> float:
     ba = b - a

tests/test_06_math/test_613_is_point_in_polygon_2d.py

@@ -1,7 +1,14 @@
 # Copyright (c) 2020-2024, Manfred Moitzi
 # License: MIT License
 import pytest
-from ezdxf.math import is_point_in_polygon_2d, Vec2, is_convex_polygon_2d
+from ezdxf.math import Vec2, is_convex_polygon_2d
+
+from ezdxf.math._construct import is_point_in_polygon_2d
+is_point_in_polygon_cy = is_point_in_polygon_2d
+try:
+    from ezdxf.acc.construct import is_point_in_polygon_2d as is_point_in_polygon_cy
+except ImportError:
+    pass
 
 
 def test_inside_horizontal_box():
@@ -157,5 +164,19 @@ def test_is_outside_cw_polygon(point: Vec2):
     assert is_point_in_polygon_2d(point, SHAPE_C_CW) == -1
 
 
+@pytest.mark.skipif(
+    is_point_in_polygon_cy is is_point_in_polygon_2d,
+    reason="no Cython implementation available",
+)
+class TestCytonImplementation:
+    @pytest.mark.parametrize("point", POINTS_INSIDE)
+    def test_is_inside_ccw_polygon(self, point: Vec2):
+        assert is_point_in_polygon_cy(point, SHAPE_C_CCW) >= 0
+
+    @pytest.mark.parametrize("point", POINTS_OUTSIDE)
+    def test_is_outside_ccw_polygon(self, point: Vec2):
+        assert is_point_in_polygon_cy(point, SHAPE_C_CCW) == -1
+
+
 if __name__ == "__main__":
     pytest.main([__file__])