ENH: Improve convex_hull execution speed.

skimage/morphology/_pnpoly.h

@@ -34,9 +34,10 @@
 extern "C" {
 #endif
 
-int pnpoly(int nr_verts, double *xp, double *yp, double x, double y)
+unsigned char pnpoly(int nr_verts, double *xp, double *yp, double x, double y)
 {
-    int i,j, c=0;
+    int i, j;
+    unsigned char c = 0;
     for (i = 0, j = nr_verts-1; i < nr_verts; j = i++) {
         if ((((yp[i]<=y) && (y<yp[j])) ||
              ((yp[j]<=y) && (y<yp[i]))) &&

skimage/morphology/_pnpoly.pyx

@@ -4,11 +4,56 @@ cimport numpy as np
 import numpy as np
 
 cdef extern from "_pnpoly.h":
+     int pnpoly(int nr_verts, double *xp, double *yp,
+                double x, double y)
+
      void npnpoly(int nr_verts, double *xp, double *yp,
                   int nr_points, double *x, double *y,
                   unsigned char *result)
 
 
+def grid_points_inside_poly(shape, verts):
+    """Test whether points on a specified grid are inside a polygon.
+
+    For each ``(r, c)`` coordinate on a grid, i.e. ``(0, 0)``, ``(0, 1)`` etc.,
+    test whether that point lies inside a polygon.
+
+    Parameters
+    ----------
+    shape : tuple (M, N)
+        Shape of the grid.
+    verts : (V, 2) array
+        Specify the V vertices of the polygon, sorted either clockwise
+        or anti-clockwise.  The first point may (but does not need to be)
+        duplicated.
+
+    Returns
+    -------
+    mask : (M, N) ndarray of bool
+        True where the grid falls inside the polygon.
+
+    """
+    cdef np.ndarray[np.double_t, ndim=1, mode="c"] vx, vy
+    verts = np.asarray(verts)
+
+    vx = verts[:, 0].astype(np.double)
+    vy = verts[:, 1].astype(np.double)
+    cdef int V = vx.shape[0]
+
+    cdef int M = shape[0]
+    cdef int N = shape[1]
+    cdef int m, n
+
+    cdef np.ndarray[dtype=np.uint8_t, ndim=2, mode="c"] out = \
+         np.zeros((M, N), dtype=np.uint8)
+
+    for m in range(M):
+        for n in range(N):
+            out[m, n] = pnpoly(V, <double*>vx.data, <double*>vy.data, m, n)
+
+    return out.view(bool)
+    
+
 def points_inside_poly(points, verts):
      """Test whether points lie inside a polygon.
 
@@ -18,6 +63,7 @@ def points_inside_poly(points, verts):
          Input points, ``(x, y)``.
      verts : (M, 2) array
          Vertices of the polygon, sorted either clockwise or anti-clockwise.
+         The first point may (but does not need to be) duplicated.
 
      Returns
      -------

skimage/morphology/convex_hull.py

@@ -2,7 +2,7 @@ __all__ = ['convex_hull']
 
 import numpy as np
 from scipy.spatial import Delaunay
-from ._pnpoly import points_inside_poly
+from ._pnpoly import points_inside_poly, grid_points_inside_poly
 
 def convex_hull(image):
     """Compute the convex hull of a binary image.
@@ -49,8 +49,6 @@ def convex_hull(image):
 
     # For each pixel coordinate, check whether that pixel
     # lies inside the convex hull
-    xy = np.dstack(np.mgrid[:image.shape[0], :image.shape[1]]).reshape(-1, 2)
-    mask = points_inside_poly(xy, v)
-    mask = mask.reshape(image.shape[:2])
+    mask = grid_points_inside_poly(image.shape[:2], v)
     
     return mask

skimage/morphology/tests/test_pnpoly.py

@@ -1,8 +1,10 @@
 import numpy as np
+from numpy.testing import assert_array_equal
 
-from skimage.morphology._pnpoly import points_inside_poly
+from skimage.morphology._pnpoly import points_inside_poly, \
+                                       grid_points_inside_poly
 
-class test_poly():
+class test_npnpoly():
     def test_square(self):
         v = np.array([[0, 0],
                       [0, 1],
@@ -22,5 +24,15 @@ class test_poly():
     def test_type(self):
         assert(points_inside_poly([[0, 0]], [[0, 0]]).dtype == np.bool)
 
+def test_grid_points_inside_poly():
+    v = np.array([[0, 0],
+                  [5, 0],
+                  [5, 5]])
+
+    expected = np.tril(np.ones((5, 5), dtype=bool))
+    
+    assert_array_equal(grid_points_inside_poly((5, 5), v),
+                       expected)
+
 if __name__ == "__main__":
     np.testing.run_module_suite()