Move points_in_polygon to Cython

This function was part of the agg C-extension. It has no dependence on anything
in AGG, so it can be moved into its own extension.

enable/_cython_speedups.pyx

@@ -0,0 +1,83 @@
+from libc.stdint cimport uint8_t
+import numpy as np
+cimport _hit_test
+
+
+def points_in_polygon(pts, poly_pts, use_winding=False):
+    """Test whether point pairs in pts are within the polygon, poly_pts.
+
+    Parameters
+    ----------
+    pts
+        an Nx2 array of x,y point pairs (floating point).  Each point is tested
+        to determine whether it falls within the polygon defined by `poly_pts`.
+    poly_pts
+        an Mx2 array of x,y point pairs (floating point) that define the
+        boundaries of a polygon.  The last point is considered to be connected
+        to the first point.
+    return
+        a 1D array of integers.  1 is returned if the corresponding x,y pair
+        in `pts` falls within `poly_pts`.  0 is returned otherwise.
+
+    This algorithm works for complex polygons.
+
+    Note: If the test point is on the border of the polygon, this
+    algorithm will deliver unpredictable results; i.e. the result
+    may be "inside" or "outside" depending on arbitrary factors
+    such as how the polygon is oriented with respect to the
+    coordinate system.
+
+    Adapted from: http://www.alienryderflex.com/polygon/
+
+    Example::
+
+        >>> import numpy as np
+        >>> from enable import points_in_polygon
+        >>> poly = np.array(((0.0, 0.0),
+                            (10.0, 0.0),
+                            (10.0, 10.0),
+                            (0.0, 10.0)))
+        >>> pts = np.array(((-1.0, -1.0),
+                           (5.0, 5.0),
+                           (15.0, 15.0)))
+        >>> results = points_in_polygon(pts, poly)
+        [0 1 0]
+
+    """
+
+    # Check the shape of pts and transpose if necessary.
+    pts = np.asarray(pts, dtype=np.float64)
+    if pts.ndim == 1:
+        pts = np.reshape(pts, (1,) + np.shape(pts))
+    if np.shape(pts)[1] != 2:
+        if np.shape(pts)[0] == 2:
+            pts = np.transpose(pts)
+        else:
+            raise ValueError('pts must be an Nx2 or 2xN array')
+
+    # Check the shape of poly_pts and transpose if necessary
+    poly_pts = np.asarray(poly_pts, dtype=np.float64)
+    if poly_pts.ndim == 1:
+        poly_pts = np.reshape(poly_pts, (1,) + np.shape(poly_pts))
+    if np.shape(poly_pts)[1] != 2:
+        if np.shape(poly_pts)[0] == 2:
+            poly_pts = np.transpose(poly_pts)
+        else:
+            raise ValueError('poly_pts must be an Nx2 or 2xN array')
+
+    cdef double[:, ::1] pts_view = pts
+    cdef double[:, ::1] poly_pts_view = poly_pts
+    cdef uint8_t[::1] results = np.zeros(len(pts), dtype=np.uint8)
+
+    if use_winding:
+        _hit_test.points_in_polygon_winding(&pts_view[0][0], pts_view.shape[0],
+                                            &poly_pts_view[0][0],
+                                            poly_pts_view.shape[0],
+                                            &results[0], results.shape[0])
+    else:
+        _hit_test.points_in_polygon(&pts_view[0][0], pts_view.shape[0],
+                                    &poly_pts_view[0][0],
+                                    poly_pts_view.shape[0],
+                                    &results[0], results.shape[0])
+
+    return results.base.astype(np.bool_)

kiva/agg/src/kiva_hit_test.cpp → enable/_hit_test.cpp

@@ -1,6 +1,6 @@
-#include "kiva_hit_test.h"
+#include "_hit_test.h"
 
-namespace kiva
+namespace enable
 {
     // Adapted from: http://www.alienryderflex.com/polygon/
     //
@@ -19,12 +19,12 @@ namespace kiva
     // coordinate system.
 
     inline bool toggle_odd_node(double x, double y,
-                            double p1x, double p1y,
-                            double p2x, double p2y)
+                                double p1x, double p1y,
+                                double p2x, double p2y)
     {
         bool toggle = false;
-        if(   ((p1y<y) && (p2y>=y))
-           || ((p2y<y) && (p1y>=y)) )
+        if ( ((p1y<y) && (p2y>=y))
+          || ((p2y<y) && (p1y>=y)) )
         {
             if (p1x + (y-p1y)/(p2y-p1y) * (p2x-p1x) < x)
             {
@@ -63,12 +63,12 @@ namespace kiva
     }
 
     void points_in_polygon(double* pts, int Npts,
-                          double* poly_pts, int Npoly_pts,
-                          unsigned char* results, int Nresults)
+                           double* poly_pts, int Npoly_pts,
+                           unsigned char* results, int Nresults)
     {
         // Nresults and Npts should match.
 
-        for(int i=0; i < Npts; i++)
+        for (int i=0; i < Npts; i++)
         {
             int ii = i*2;
             double x = pts[ii];
@@ -77,12 +77,16 @@ namespace kiva
         }
     }
 
-    inline double is_left(double x, double y, double x1, double y1, double x2, double y2)
+    inline double is_left(double x, double y,
+                          double x1, double y1,
+                          double x2, double y2)
     {
         return ( (x2 - x1) * (y - y1) - (x - x1) * (y2 - y1) );
     }
 
-    inline int winding_increment(double x, double y, double x1, double y1, double x2, double y2)
+    inline int winding_increment(double x, double y,
+                                 double x1, double y1,
+                                 double x2, double y2)
     {
         if (y1 <= y)
         {
@@ -108,7 +112,8 @@ namespace kiva
     }
 
 
-    bool point_in_polygon_winding(double x, double y, double* poly_pts, int Npoly_pts)
+    bool point_in_polygon_winding(double x, double y,
+                                  double* poly_pts, int Npoly_pts)
     {
         int winding_number = 0;
         double p1_x, p1_y, p2_x, p2_y;
@@ -120,7 +125,7 @@ namespace kiva
             p2_x = poly_pts[ii+2];
             p2_y = poly_pts[ii+3];
 
-            winding_number += winding_increment(x,y,p1_x,p1_y,p2_x,p2_y);
+            winding_number += winding_increment(x, y, p1_x, p1_y, p2_x, p2_y);
         }
 
         // Last point wraps to the beginning.
@@ -129,18 +134,18 @@ namespace kiva
         p2_x = poly_pts[0];
         p2_y = poly_pts[1];
 
-        winding_number += winding_increment(x,y,p1_x,p1_y,p2_x,p2_y);
+        winding_number += winding_increment(x, y, p1_x, p1_y, p2_x, p2_y);
 
         return winding_number != 0;
     }
 
     void points_in_polygon_winding(double* pts, int Npts,
-                          double* poly_pts, int Npoly_pts,
-                          unsigned char* results, int Nresults)
+                                   double* poly_pts, int Npoly_pts,
+                                   unsigned char* results, int Nresults)
     {
         // Nresults and Npts should match.
 
-        for(int i=0; i < Npts; i++)
+        for (int i=0; i < Npts; i++)
         {
             int ii = i*2;
             double x = pts[ii];

kiva/agg/src/kiva_hit_test.h → enable/_hit_test.h

@@ -1,7 +1,7 @@
-#ifndef KIVA_HIT_TEST_H
-#define KIVA_HIT_TEST_H
+#ifndef ENABLE_HIT_TEST_H
+#define ENABLE_HIT_TEST_H
 
-namespace kiva
+namespace enable
 {
     bool point_in_polygon(double x, double y, 
                           double* poly_pts, int Npoly_pts);
@@ -13,6 +13,6 @@ namespace kiva
     void points_in_polygon_winding(double* pts, int Npts, 
                                    double* poly_pts, int Npoly_pts,
                                    unsigned char* results, int Nresults);
-	    		  
+
 }
 #endif

enable/_hit_test.pxd

@@ -0,0 +1,10 @@
+
+
+cdef extern from "_hit_test.h" namespace "enable":
+    void points_in_polygon(double* pts, int Npts,
+                          double* poly_pts, int Npoly_pts,
+                          unsigned char* results, int Nresults)
+
+    void points_in_polygon_winding(double* pts, int Npts,
+                                   double* poly_pts, int Npoly_pts,
+                                   unsigned char* results, int Nresults)

enable/api.py

@@ -63,4 +63,6 @@
 from .viewport import Viewport
 from .window import Window
 
+from ._cython_speedups import points_in_polygon
+
 from .primitives.api import Annotater, Box, Line, Polygon

enable/primitives/polygon.py

@@ -6,13 +6,12 @@
 
 # Enthought library imports.
 from kiva.constants import EOF_FILL_STROKE, FILL, FILL_STROKE
-from kiva.agg import points_in_polygon
-from traits.api import Any, Event, Float, HasTraits, Instance, List, \
-                             Property, Trait, Tuple
+from traits.api import (Any, Event, Float, HasTraits, Instance, List,
+                        Property, Trait, Tuple)
 from traitsui.api import Group, View
 
 # Local imports.
-from enable.api import border_size_trait, Component
+from enable.api import border_size_trait, points_in_polygon, Component
 from enable.colors import ColorTrait
 
 

kiva/agg/tests/points_in_polygon_test_case.py → enable/tests/points_in_polygon_test.py

@@ -2,7 +2,8 @@
 
 from numpy import array, allclose
 
-from kiva import agg
+from enable.api import points_in_polygon
+
 
 class TestPointsInPolygon(unittest.TestCase):
 
@@ -11,8 +12,8 @@ def test_simple_points_in_polygon(self):
         polygon = array(((0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0)))
         points = array(((-1.0, -1.0), (5.0, 5.0), (15.0, 15.0)))
 
-        result = agg.points_in_polygon(points, polygon)
-        self.assertTrue(allclose(array([0,1,0]), result))
+        result = points_in_polygon(points, polygon)
+        self.assertTrue(allclose(array([0, 1, 0]), result))
 
         return
 
@@ -21,38 +22,37 @@ def test_asymmetric_points_in_polygon(self):
         polygon = array(((0.0, 0.0), (20.0, 0.0), (20.0, 10.0), (0.0, 10.0)))
         points = array(((5.0, 5.0), (10.0, 5.0), (15.0, 5.0)))
 
-        result = agg.points_in_polygon(points, polygon)
-        self.assertTrue(allclose(array([1,1,1]), result))
+        result = points_in_polygon(points, polygon)
+        self.assertTrue(allclose(array([1, 1, 1]), result))
 
         return
 
-
     def test_rectangle(self):
 
-        vertices = array(((0,0), (0,10), (10,10), (10,0)))
+        vertices = array(((0, 0), (0, 10), (10, 10), (10, 0)))
 
         # Try the lower left.
-        trial = array(((0,0),))
-        oe_result = agg.points_in_polygon(trial, vertices)
-        w_result = agg.points_in_polygon(trial, vertices, True)
+        trial = array(((0, 0),))
+        oe_result = points_in_polygon(trial, vertices)
+        w_result = points_in_polygon(trial, vertices, True)
         self.assertEqual(0, oe_result[0],
                          "Lower left corner not in polygon. OEF")
         self.assertEqual(1, w_result[0],
                          "Lower left corner not in polygon. Winding")
 
         # Try the center.
-        trial = array(((5,5),))
-        oe_result = agg.points_in_polygon(trial, vertices)
-        w_result = agg.points_in_polygon(trial, vertices, True)
+        trial = array(((5, 5),))
+        oe_result = points_in_polygon(trial, vertices)
+        w_result = points_in_polygon(trial, vertices, True)
         self.assertEqual(1, oe_result[0],
                          "Center not in polygon. OEF")
         self.assertEqual(1, w_result[0],
                          "Center not in polygon. Winding")
 
         # Try the center.
-        trial = array(((10,10),))
-        oe_result = agg.points_in_polygon(trial, vertices)
-        w_result = agg.points_in_polygon(trial, vertices, True)
+        trial = array(((10, 10),))
+        oe_result = points_in_polygon(trial, vertices)
+        w_result = points_in_polygon(trial, vertices, True)
         self.assertEqual(1, oe_result[0],
                          "Top-right in polygon. OEF")
         self.assertEqual(0, w_result[0],
@@ -78,7 +78,7 @@ def test_center_removed(self):
         # The inner square with containing the edge (5,6) is outside in both
         # cases.
 
-        vertices = array(((0,0),
+        vertices = array(((0, 0),
                           (10, 0),
                           (10, 8),
                           (2, 8),
@@ -89,24 +89,23 @@ def test_center_removed(self):
                           (5, 10),
                           (0, 10)))
 
-        trial = array(((3,7),))
-        oe_result = agg.points_in_polygon(trial, vertices)
-        w_result = agg.points_in_polygon(trial, vertices, True)
+        trial = array(((3, 7),))
+        oe_result = points_in_polygon(trial, vertices)
+        w_result = points_in_polygon(trial, vertices, True)
         self.assertEqual(0, oe_result[0],
                          "Interior polygon inside: odd-even")
         self.assertEqual(1, w_result[0],
                          "Interior polygon outside: winding")
 
-        trial = array(((6,5),))
-        oe_result = agg.points_in_polygon(trial, vertices)
-        w_result = agg.points_in_polygon(trial, vertices, True)
+        trial = array(((6, 5),))
+        oe_result = points_in_polygon(trial, vertices)
+        w_result = points_in_polygon(trial, vertices, True)
         self.assertEqual(0, oe_result[0],
                          "Interior polygon inside: odd-even")
         self.assertEqual(0, w_result[0],
                          "Interior polygon inside: winding")
 
 
 if __name__ == "__main__":
-          unittest.main()
+    unittest.main()
 
-#### EOF ######################################################################

kiva/agg/agg.i

@@ -1,10 +1,10 @@
 /* -*- c++ -*- */
 /* File : example.i */
 %module agg
-
-#if (SWIG_VERSION > 0x010322)
-%feature("compactdefaultargs");
-#endif // (SWIG_VERSION > 0x010322)
+
+#if (SWIG_VERSION > 0x010322)
+%feature("compactdefaultargs");
+#endif // (SWIG_VERSION > 0x010322)
 
 %include "constants.i"
 %include "rgba.i"
@@ -13,5 +13,4 @@
 %include "compiled_path.i"
 //%include "image.i"
 %include "graphics_context.i"
-%include "hit_test.i"