ENH: Only perform convex hull on select candidate pixels. Gets us a factor 30 speedup.

skimage/morphology/_convex_hull.pyx

@@ -0,0 +1,59 @@
+# -*- python -*-
+
+cimport numpy as np
+import numpy as np
+
+def possible_hull(np.ndarray[dtype=np.uint8_t, ndim=2, mode="c"] img):
+    """Return positions of pixels that possibly belong to the convex hull.
+
+    Parameters
+    ----------
+    img : ndarray of bool
+        Binary input image.
+
+    Returns
+    -------
+    coords : ndarray (N, 2)
+       The ``(row, column)`` coordinates of all pixels that possibly belong to
+       the convex hull.
+
+    """
+    cdef int i, j, k
+    cdef unsigned int M, N
+    
+    M = img.shape[0]
+    N = img.shape[1]
+
+    # Output: M storage slots for left boundary pixels
+    #         N storage slots for top boundary pixels
+    #         M storage slots for right boundary pixels
+    #         N storage slots for bottom boundary pixels
+    cdef np.ndarray[dtype=np.int_t, ndim=2] nonzero = \
+         np.ones((2 * (M + N), 2), dtype=np.int)
+    nonzero *= -1 
+
+    k = 0
+    for i in range(M):
+        for j in range(N):
+            if img[i, j] != 0:
+                # Left check
+                if nonzero[i, 1] == -1:
+                    nonzero[i, 0] = i
+                    nonzero[i, 1] = j
+
+                # Right check
+                elif nonzero[M + N + i, 1] < j:
+                    nonzero[M + N + i, 0] = i
+                    nonzero[M + N + i, 1] = j
+
+                # Top check
+                if nonzero[M + j, 1] == -1:
+                    nonzero[M + j, 0] = i
+                    nonzero[M + j, 1] = j
+
+                # Bottom check
+                elif nonzero[2 * M + N + j, 0] < i:
+                    nonzero[2 * M + N + j, 0] = i
+                    nonzero[2 * M + N + j, 1] = j
+    
+    return nonzero[nonzero[:, 0] != -1]

skimage/morphology/convex_hull.py

@@ -3,6 +3,7 @@ __all__ = ['convex_hull']
 import numpy as np
 from scipy.spatial import Delaunay
 from ._pnpoly import points_inside_poly, grid_points_inside_poly
+from ._convex_hull import possible_hull
 
 def convex_hull(image):
     """Compute the convex hull of a binary image.
@@ -25,11 +26,16 @@ def convex_hull(image):
     .. [1] http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/
 
     """
-    image = image.astype(bool)
-    r, c = np.nonzero(image)
-    coords = np.vstack((r, c)).T
-    N = len(coords)
 
+    image = image.astype(bool)
+
+    # Here we do an optimisation by choosing only pixels that are
+    # the starting or ending pixel of a row or column.  This vastly
+    # limits the number of coordinates to examine for the virtual
+    # hull.
+    coords = possible_hull(image.astype(np.uint8))
+    N = len(coords)
+    
     # Add a vertex for the middle of each pixel edge
     coords_corners = np.empty((N * 4, 2))
     for i, (x_offset, y_offset) in enumerate(zip((0, 0, -0.5, 0.5),

skimage/morphology/setup.py

@@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None):
     cython(['_watershed.pyx'], working_path=base_path)
     cython(['_skeletonize.pyx'], working_path=base_path)
     cython(['_pnpoly.pyx'], working_path=base_path)
+    cython(['_convex_hull.pyx'], working_path=base_path)
 
     config.add_extension('ccomp', sources=['ccomp.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -27,7 +28,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_pnpoly', sources=['_pnpoly.c'],
                          include_dirs=[get_numpy_include_dirs()])
-
+    config.add_extension('_convex_hull', sources=['_convex_hull.c'],
+                         include_dirs=[get_numpy_include_dirs()])
 
     return config
 

skimage/morphology/tests/test_convex_hull.py

@@ -1,6 +1,7 @@
 import numpy as np
 from numpy.testing import assert_array_equal
 from skimage.morphology import convex_hull
+from skimage.morphology._convex_hull import possible_hull
 
 def test_basic():
     image = np.array(
@@ -21,5 +22,38 @@ def test_basic():
 
     assert_array_equal(convex_hull(image), expected)
 
+
+def test_possible_hull():
+    image = np.array(
+        [[0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 1, 0, 0, 0, 0],
+         [0, 0, 0, 1, 0, 1, 0, 0, 0],
+         [0, 0, 1, 1, 1, 1, 1, 0, 0],
+         [0, 1, 1, 1, 1, 1, 1, 1, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8)
+
+    expected = np.array([[1, 4],
+                         [2, 3],
+                         [3, 2],
+                         [4, 1],
+                         [4, 1],
+                         [3, 2],
+                         [2, 3],
+                         [1, 4],
+                         [2, 5],
+                         [3, 6],
+                         [4, 7],
+                         [2, 5],
+                         [3, 6],
+                         [4, 7],
+                         [4, 2],
+                         [4, 3],
+                         [4, 4],
+                         [4, 5],
+                         [4, 6]])
+
+    ph = possible_hull(image)
+    assert_array_equal(possible_hull(image), expected)
+
 if __name__ == "__main__":
     np.testing.run_module_suite()