Add convex hull and point in polygon functions

2 files changed, 82 insertions, 0 deletions

diff --git a/gerbonara/tests/test_utils.py b/gerbonara/tests/test_utils.py
index 1791afe..bc3fedb 100644
--- a/gerbonara/tests/test_utils.py
+++ b/gerbonara/tests/test_utils.py
@@ -17,8 +17,14 @@
 # limitations under the License.
 #
 
+import math
+import random
+
 import pytest
+
 from ..cam import FileSettings
+from ..utils import convex_hull, point_in_polygon, setup_svg, Tag
+from .utils import *
 
 
 def test_zero_suppression():
@@ -103,3 +109,25 @@ def test_write_format_validation():
             settings = FileSettings(number_format=fmt)
             settings.write_gerber_value(69.0)
 
+def test_convex_hull_and_point_in_polygon(tmpfile):
+    svg = tmpfile('Visualization', '.svg')
+    st = random.Random(0)
+    for _ in range(50):
+        for n in [*range(1, 10), 12, 15, 20, 30, 50, 300, 1000, 5000]:
+            w = math.sqrt(n) * 10
+            rd = lambda: round(st.random() * w)
+            rp = lambda: (rd(), rd())
+            points = {rp() for _ in range(n)}
+            hull_l = convex_hull(points)
+            hull = set(hull_l)
+
+            tags = [Tag('circle', cx=x, cy=y, r=0.2, fill=('red' if (x, y) in hull else 'black')) for x, y in points]
+            for (x0, y0), (x1, y1) in zip([hull_l[-1], *hull_l[:-1]], hull_l):
+                tags.append(Tag('path', d=f'M {x0},{y0} L {x1},{y1}', stroke_width='0.1', stroke='red', fill='none'))
+            svg.write_text(str(setup_svg(tags, bounds=((0, 0), (w, w)), margin=1)))
+
+            # all hull corners must be in the set of original points
+            assert not (hull-points)
+            for p in points-hull:
+                assert point_in_polygon(p, hull_l)
+
diff --git a/gerbonara/utils.py b/gerbonara/utils.py
index 6b2d5c1..53f6398 100644
--- a/gerbonara/utils.py
+++ b/gerbonara/utils.py
@@ -28,6 +28,7 @@ This module provides utility functions for working with Gerber and Excellon file
 import os
 import re
 import textwrap
+from functools import reduce
 from enum import Enum
 import math
 
@@ -396,6 +397,33 @@ def arc_bounds(x1, y1, x2, y2, cx, cy, clockwise):
     return (min_x+cx, min_y+cy), (max_x+cx, max_y+cy)
 
 
+def convex_hull(points):
+    '''
+    Returns points on convex hull in CCW order according to Graham's scan algorithm. 
+    By Tom Switzer <thomas.switzer@gmail.com>.
+    '''
+    # https://gist.github.com/arthur-e/5cf52962341310f438e96c1f3c3398b8
+    TURN_LEFT, TURN_RIGHT, TURN_NONE = (1, -1, 0)
+
+    def cmp(a, b):
+        return (a > b) - (a < b)
+
+    def turn(p, q, r):
+        return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0)
+
+    def keep_left(hull, r):
+        while len(hull) > 1 and turn(hull[-2], hull[-1], r) != TURN_LEFT:
+            hull.pop()
+        if not len(hull) or hull[-1] != r:
+            hull.append(r)
+        return hull
+
+    points = sorted(points)
+    l = reduce(keep_left, points, [])
+    u = reduce(keep_left, reversed(points), [])
+    return l.extend(u[i] for i in range(1, len(u) - 1)) or l
+
+
 def point_line_distance(l1, l2, p):
     """ Calculate distance between infinite line through l1 and l2, and point p. """
     # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
@@ -471,3 +499,29 @@ def setup_svg(tags, bounds, margin=0, arg_unit=MM, svg_unit=MM, pagecolor='white
             **namespaces,
             root=True)
 
+
+def point_in_polygon(point, poly):
+    # https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
+    # https://wrfranklin.org/Research/Short_Notes/pnpoly.html
+
+    if not poly:
+        return False
+
+    res = False
+    tx, ty = point
+    xp, yp = poly[-1]
+    for x, y in poly:
+        if yp == ty == y and ((x > tx) != (xp > tx)): # test point on horizontal segment
+            return True
+        if xp == tx == x and ((y > ty) != (yp > ty)): # test point on vertical segment
+            return True
+        if ((y > ty) != (yp > ty)):
+            tmp = ((xp-x) * (ty-y) / (yp-y) + x)
+            if tx == tmp: # test point on diagonal segment
+                return True
+            elif tx < tmp:
+                res = not res
+        xp, yp = x, y
+
+    return res
+