draft: add c extension numpy

.github/workflows/pytest.yml

@@ -25,6 +25,7 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
+        python setup.py build_ext --inplace
         python -m pip install .
         python -m pip install pytest
     - name: Test with pytest

boundvor/geometry.py

@@ -1,78 +1,29 @@
 import numpy as np
+import boundvor.geometry_c as c_api
 
 def point_in_polygon(point, polygon, check_boundary=False):
 
     if len(polygon) < 3:
-        raise ValueError("Polygon must have at least three points")
+         raise ValueError("Polygon must have at least three points")
 
-    x, y = point
-    inside = False
-    n = len(polygon)
-    p1x, p1y = polygon[0]
-    for i in range(n + 1):
-        p2x, p2y = polygon[i % n]
-        if y > min(p1y, p2y):
-            if y <= max(p1y, p2y):
-                if x <= max(p1x, p2x):
-                    if p1y != p2y:
-                        xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
-                    if p1x == p2x or x <= xinters:
-                        inside = not inside
-        p1x, p1y = p2x, p2y
+    point = np.asarray(point, dtype=np.float64)
+    polygon = np.asarray(polygon, dtype=np.float64)
 
-    if check_boundary:
-        for i in range(n):
-            p1x, p1y = polygon[i]
-            p2x, p2y = polygon[(i + 1) % n]
-            if (min(p1x, p2x) <= x <= max(p1x, p2x)) and (min(p1y, p2y) <= y <= max(p1y, p2y)):
-                if (p2x - p1x) * (y - p1y) == (p2y - p1y) * (x - p1x):
-                    return True
+    return c_api.point_in_polygon(point, polygon, check_boundary)
 
-    return inside
 
-def line_intersect(p1, p2, q1, q2):
-    def ccw(A, B, C):
-        return (C[1] - A[1]) * (B[0] - A[0]) > (B[1] - A[1]) * (C[0] - A[0])
-    return ccw(p1, q1, q2) != ccw(p2, q1, q2) and ccw(p1, p2, q1) != ccw(p1, p2, q2)
+def _sort_counter_clockwise(polygon):
+    (center_x, center_y) = np.mean(polygon, axis=0)
+    key_fun = lambda x: np.arctan2(x[1] - center_y, x[0] - center_x)
+    return sorted(polygon, key=key_fun)
 
-def polygon_intersection(poly1, poly2):
-    def clip(subject_polygon, clip_polygon):
-        def inside(p, edge):
-            return (edge[1][0] - edge[0][0]) * (p[1] - edge[0][1]) >= (edge[1][1] - edge[0][1]) * (p[0] - edge[0][0])
 
-        def compute_intersection(s, e, edge):
-            dc = [edge[0][0] - edge[1][0], edge[0][1] - edge[1][1]]
-            dp = [s[0] - e[0], s[1] - e[1]]
-            n1 = edge[0][0] * edge[1][1] - edge[0][1] * edge[1][0]
-            n2 = s[0] * e[1] - s[1] * e[0]
-            n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0])
-            return [(n1 * dp[0] - n2 * dc[0]) * n3, (n1 * dp[1] - n2 * dc[1]) * n3]
+def polygon_intersection(subject_polygon, clipping_polygon):
 
-        output_list = subject_polygon
-        clip_polygon.append(clip_polygon[0])  # Ensure the clip polygon is closed
-        for clip_edge in zip(clip_polygon[:-1], clip_polygon[1:]):
-            input_list = output_list
-            output_list = []
-            if len(input_list) == 0:
-                break
-            s = input_list[-1]
-            for e in input_list:
-                if inside(e, clip_edge):
-                    if not inside(s, clip_edge):
-                        output_list.append(compute_intersection(s, e, clip_edge))
-                    output_list.append(e)
-                elif inside(s, clip_edge):
-                    output_list.append(compute_intersection(s, e, clip_edge))
-                s = e
-        return output_list
-
-    def sort_polygon(poly):
-        center = np.mean(poly, axis=0)
-        return sorted(poly, key=lambda x: np.arctan2(x[1] - center[1], x[0] - center[0]))
-
-    poly1 = sort_polygon(poly1)
-    poly2 = sort_polygon(poly2)
-
-    clipped = clip(poly1, poly2)
-    return clipped
+    # Sort counter-clockwise
+    subject_polygon = _sort_counter_clockwise(list(subject_polygon))
+    clipping_polygon = _sort_counter_clockwise(list(clipping_polygon))
 
+    clipped_polygon = c_api.polygon_intersection(subject_polygon, clipping_polygon)
+    clipped_polygon = np.asarray(clipped_polygon, dtype=np.float64)
+    return clipped_polygon

pyproject.toml

@@ -1,5 +1,5 @@
 [build-system]
-requires = ["setuptools", "wheel"]
+requires = ["setuptools", "wheel", "numpy"]
 build-backend = "setuptools.build_meta"
 
 [project]

setup.py

@@ -0,0 +1,15 @@
+from setuptools import setup, Extension
+import numpy
+
+def get_numpy_include_dirs():
+    return [numpy.get_include()]
+
+geometry_c_extension = Extension(
+    "boundvor.geometry_c",
+    sources=["src/geometry.c"],
+    include_dirs=get_numpy_include_dirs()
+)
+
+setup(
+    ext_modules=[geometry_c_extension],
+)

src/geometry.c

@@ -0,0 +1,224 @@
+#define PY_SSIZE_T_CLEAN
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#include <Python.h>
+#include <numpy/arrayobject.h>
+#include <math.h>
+
+static int point_in_polygon(double x, double y, double *polygon, int n, int check_boundary) {
+    int i, j;
+    int inside = 0;
+
+    // Check if the point is inside the polygon
+    for (i = 0, j = n - 1; i < n; j = i++) {
+        double xi = polygon[2 * i];
+        double yi = polygon[2 * i + 1];
+        double xj = polygon[2 * j];
+        double yj = polygon[2 * j + 1];
+
+        int intersect = ((yi > y) != (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi);
+        if (intersect) {
+            inside = !inside;
+        }
+    }
+
+    // Check if the point is on the boundary of the polygon
+    if (check_boundary) {
+        for (i = 0, j = n - 1; i < n; j = i++) {
+            double xi = polygon[2 * i];
+            double yi = polygon[2 * i + 1];
+            double xj = polygon[2 * j];
+            double yj = polygon[2 * j + 1];
+            if ((x >= fmin(xi, xj)) && (x <= fmax(xi, xj)) && (y >= fmin(yi, yj)) && (y <= fmax(yi, yj))) {
+                if ((xj - xi) * (y - yi) == (yj - yi) * (x - xi)) {
+                    return 1;  // Point is on the boundary
+                }
+            }
+        }
+    }
+
+    return inside;
+}
+
+// Helper function to determine if a point is inside a clipping edge
+static int inside(double *p, double *edge_start, double *edge_end) {
+    return (edge_end[0] - edge_start[0]) * (p[1] - edge_start[1]) >= (edge_end[1] - edge_start[1]) * (p[0] - edge_start[0]);
+}
+
+static void compute_intersection(double *s, double *e, double *edge_start, double *edge_end, double *result) {
+    double dc[2] = {edge_start[0] - edge_end[0], edge_start[1] - edge_end[1]};
+    double dp[2] = {s[0] - e[0], s[1] - e[1]};
+    double n1 = edge_start[0] * edge_end[1] - edge_start[1] * edge_end[0];
+    double n2 = s[0] * e[1] - s[1] * e[0];
+    double det = dc[0] * dp[1] - dc[1] * dp[0];
+
+    if (det == 0) {
+        result[0] = result[1] = NAN;
+        return;
+    }
+
+    result[0] = (n1 * dp[0] - n2 * dc[0]) / det;
+    result[1] = (n1 * dp[1] - n2 * dc[1]) / det;
+}
+
+
+// Computes the interesection between two polygons using the Sutherland-Hodgman algorithm
+static PyObject* polygon_intersection(double *subject_polygon, int subject_size, double *clip_polygon, int clip_size) {
+    PyObject *result = PyList_New(0);
+    if (!result) return NULL;
+
+    // Allocate memory for closed clip polygon (add the first point at the end)
+    double *clip_polygon_closed = malloc((clip_size + 1) * 2 * sizeof(double));
+    if (!clip_polygon_closed) {
+        Py_DECREF(result);
+        return NULL;
+    }
+    memcpy(clip_polygon_closed, clip_polygon, clip_size * 2 * sizeof(double));
+    clip_polygon_closed[2 * clip_size] = clip_polygon[0];
+    clip_polygon_closed[2 * clip_size + 1] = clip_polygon[1];
+
+    // Allocate memory for output list
+    double *output_list = malloc(subject_size * 2 * sizeof(double));
+    if (!output_list) {
+        free(clip_polygon_closed);
+        Py_DECREF(result);
+        return NULL;
+    }
+    memcpy(output_list, subject_polygon, subject_size * 2 * sizeof(double));
+    int output_size = subject_size * 2;
+
+    for (int i = 0; i < clip_size; ++i) {
+        double *clip_edge_start = clip_polygon_closed + 2 * i;
+        double *clip_edge_end = clip_polygon_closed + 2 * (i + 1);
+
+        double *input_list = malloc(output_size * sizeof(double));
+        if (!input_list) {
+            free(clip_polygon_closed);
+            free(output_list);
+            Py_DECREF(result);
+            return NULL;
+        }
+        memcpy(input_list, output_list, output_size * sizeof(double));
+        int input_size = output_size / 2;
+        output_size = 0;
+
+        double *s = input_list + 2 * (input_size - 1);
+        for (int j = 0; j < input_size; ++j) {
+            double *e = input_list + 2 * j;
+
+            double intersection[2];
+            compute_intersection(s, e, clip_edge_start, clip_edge_end, intersection);
+
+            if (inside(e, clip_edge_start, clip_edge_end)) {
+                if (!inside(s, clip_edge_start, clip_edge_end)) {
+                    output_list[output_size++] = intersection[0];
+                    output_list[output_size++] = intersection[1];
+                }
+                output_list[output_size++] = e[0];
+                output_list[output_size++] = e[1];
+            } else if (inside(s, clip_edge_start, clip_edge_end)) {
+                output_list[output_size++] = intersection[0];
+                output_list[output_size++] = intersection[1];
+            }
+            s = e;
+        }
+        free(input_list);
+    }
+
+    free(clip_polygon_closed);
+
+    for (int i = 0; i < output_size; i += 2) {
+        PyObject *point = Py_BuildValue("(dd)", output_list[i], output_list[i + 1]);
+        if (!point) {
+            // If Py_BuildValue fails, clean up and return NULL
+            free(output_list);
+            Py_DECREF(result);
+            return NULL;
+        }
+        PyList_Append(result, point);
+        Py_DECREF(point);
+    }
+
+    free(output_list);
+    return result;
+}
+
+static PyObject* py_point_in_polygon(PyObject* self, PyObject* args) {
+    PyObject *point_obj, *polygon_obj;
+    int check_boundary;
+
+    if (!PyArg_ParseTuple(args, "OOi", &point_obj, &polygon_obj, &check_boundary)) {
+        return NULL;
+    }
+
+    PyArrayObject *point_array = (PyArrayObject*) PyArray_FROM_OTF(point_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    PyArrayObject *polygon_array = (PyArrayObject*) PyArray_FROM_OTF(polygon_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+
+    if (!point_array || !polygon_array) {
+        Py_XDECREF(point_array);
+        Py_XDECREF(polygon_array);
+        return NULL;
+    }
+
+    double x = *(double*) PyArray_GETPTR1(point_array, 0);
+    double y = *(double*) PyArray_GETPTR1(point_array, 1);
+    double *polygon = (double*) PyArray_DATA(polygon_array);
+    int n = PyArray_DIM(polygon_array, 0);
+
+    int result = point_in_polygon(x, y, polygon, n, check_boundary);
+
+    Py_DECREF(point_array);
+    Py_DECREF(polygon_array);
+
+    return PyBool_FromLong(result);
+}
+
+static PyObject* py_polygon_intersection(PyObject* self, PyObject* args) {
+    PyObject *poly1_obj, *poly2_obj;
+    if (!PyArg_ParseTuple(args, "OO", &poly1_obj, &poly2_obj)) {
+        return NULL;
+    }
+
+    PyArrayObject *poly1_array = (PyArrayObject*) PyArray_FROM_OTF(poly1_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    PyArrayObject *poly2_array = (PyArrayObject*) PyArray_FROM_OTF(poly2_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+
+    if (!poly1_array || !poly2_array) {
+        Py_XDECREF(poly1_array);
+        Py_XDECREF(poly2_array);
+        return NULL;
+    }
+
+    double *poly1 = (double*) PyArray_DATA(poly1_array);
+    int poly1_size = PyArray_DIM(poly1_array, 0);
+    double *poly2 = (double*) PyArray_DATA(poly2_array);
+    int poly2_size = PyArray_DIM(poly2_array, 0);
+
+    PyObject *result = polygon_intersection(poly1, poly1_size, poly2, poly2_size);
+
+    Py_DECREF(poly1_array);
+    Py_DECREF(poly2_array);
+
+    return result;
+}
+
+// Method definitions
+static PyMethodDef GeometryMethods[] = {
+    {"point_in_polygon", py_point_in_polygon, METH_VARARGS, "Check if a point is inside a polygon"},
+    {"polygon_intersection", py_polygon_intersection, METH_VARARGS, "Compute polygon intersection"},
+    {NULL, NULL, 0, NULL}
+};
+
+// Module definition
+static struct PyModuleDef geometrymodule = {
+    PyModuleDef_HEAD_INIT,
+    "geometry_c",
+    NULL,
+    -1,
+    GeometryMethods
+};
+
+// Module initialization function
+PyMODINIT_FUNC PyInit_geometry_c(void) {
+    PyObject *module = PyModule_Create(&geometrymodule);
+    import_array();
+    return module;
+}

tests/test_voronoi.py

@@ -38,8 +38,8 @@ def test_default_bounds():
 
 
 def test_furthest_site():
-    points = np.random.rand(10, 2)
-    bounding_box = np.array([[0.0, 0.0], [0.0, 1.0], [1.0, 1.0], [1.0, 0.0]])
+    points = np.array([[0.3, 0.3], [0.7, 0.7], [0.5, 0.5], [0.4, 0.2]])
+    bounding_box = [(0., 0.), (0., 1.), (1., 1.), (1., 0.)]
     voronoi = BoundedVoronoi(points, bounds=bounding_box, furthest_site=True)
 
     assert len(voronoi.regions) == len(points)