From b5ad567a686a93900e7c52e0ec61ff0715dc32e2 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Wed, 26 Oct 2011 17:42:04 -0700
Subject: [PATCH 1/7] ENH: Add binary convex hull computation.

---
 CONTRIBUTORS.txt                             |  2 +
 doc/examples/plot_convex_hull.py             | 35 ++++++++++
 skimage/morphology/__init__.py               |  1 +
 skimage/morphology/_pnpoly.h                 | 71 ++++++++++++++++++++
 skimage/morphology/_pnpoly.pyx               | 47 +++++++++++++
 skimage/morphology/convex_hull.py            | 56 +++++++++++++++
 skimage/morphology/setup.py                  |  3 +
 skimage/morphology/tests/test_convex_hull.py | 25 +++++++
 skimage/morphology/tests/test_pnpoly.py      | 26 +++++++
 9 files changed, 266 insertions(+)
 create mode 100644 doc/examples/plot_convex_hull.py
 create mode 100644 skimage/morphology/_pnpoly.h
 create mode 100644 skimage/morphology/_pnpoly.pyx
 create mode 100644 skimage/morphology/convex_hull.py
 create mode 100644 skimage/morphology/tests/test_convex_hull.py
 create mode 100644 skimage/morphology/tests/test_pnpoly.py

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 307379f8..bff86ccc 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -84,3 +84,5 @@
 - Nelle Varoquaux
   Renaming of the package to ``skimage``.
 
+- W. Randolph Franklin
+  Point in polygon test.
diff --git a/doc/examples/plot_convex_hull.py b/doc/examples/plot_convex_hull.py
new file mode 100644
index 00000000..8308bd90
--- /dev/null
+++ b/doc/examples/plot_convex_hull.py
@@ -0,0 +1,35 @@
+"""
+===========
+Convex Hull
+===========
+
+The convex hull of a binary image is the set of pixels included in the
+smallest convex polygon that surround all white pixels in the input.
+
+In this example, we show how the input pixels (white) get filled in by the
+convex hull (white and grey).
+
+A good overview of the algorithm is given on `Steve Eddin's blog
+<http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/>`__.
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage.morphology import convex_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, 0, 0, 0, 1, 0, 0],
+     [0, 1, 0, 0, 0, 0, 0, 1, 0],
+     [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=float)
+
+chull = convex_hull(image)
+image[chull] += 1.7
+image -= -1.7
+
+plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
+plt.show()
diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py
index efb92507..5a7d67d7 100644
--- a/skimage/morphology/__init__.py
+++ b/skimage/morphology/__init__.py
@@ -3,3 +3,4 @@ from selem import *
 from .ccomp import label
 from watershed import watershed, is_local_maximum
 from skeletonize import skeletonize, medial_axis
+from .convex_hull import convex_hull
diff --git a/skimage/morphology/_pnpoly.h b/skimage/morphology/_pnpoly.h
new file mode 100644
index 00000000..a477a70f
--- /dev/null
+++ b/skimage/morphology/_pnpoly.h
@@ -0,0 +1,71 @@
+/* `pnpoly` is from
+
+   http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
+
+   Copyright (c) 1970-2003, Wm. Randolph Franklin
+
+   Permission is hereby granted, free of charge, to any person
+   obtaining a copy of this software and associated documentation
+   files (the "Software"), to deal in the Software without
+   restriction, including without limitation the rights to use, copy,
+   modify, merge, publish, distribute, sublicense, and/or sell copies
+   of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+    1. Redistributions of source code must retain the above copyright
+       notice, this list of conditions and the following disclaimers.
+    2. Redistributions in binary form must reproduce the above
+       copyright notice in the documentation and/or other materials
+       provided with the distribution.
+    3. The name of W. Randolph Franklin may not be used to endorse or
+       promote products derived from this Software without specific
+       prior written permission.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+   SOFTWARE. */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int pnpoly(int nr_verts, double *xp, double *yp, double x, double y)
+{
+    int i,j, 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]))) &&
+            (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
+	    
+	    c = !c;
+    }
+    return c;
+}
+
+void npnpoly(int nr_verts, double *xp, double *yp,
+	     int nr_points, double *x, double *y,
+	     unsigned char *result)
+/*
+ * For N provided points, calculate whether they are in
+ * the polygon defined by vertices *xp, *yp.
+ *
+ * nr_verts : number of vertices
+ * *xp, *yp : x and y coordinates of vertices
+ * nr_points : number of data points provided
+ * *x, *y : data points
+ */
+{
+    unsigned char n = 0;
+    for (n = 0; n < nr_points; n++) {
+	result[n] = pnpoly(nr_verts, xp, yp, x[n], y[n]);
+    }
+}
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/skimage/morphology/_pnpoly.pyx b/skimage/morphology/_pnpoly.pyx
new file mode 100644
index 00000000..65c9e55e
--- /dev/null
+++ b/skimage/morphology/_pnpoly.pyx
@@ -0,0 +1,47 @@
+# -*- python -*-
+
+cimport numpy as np
+import numpy as np
+
+cdef extern from "_pnpoly.h":
+     void npnpoly(int nr_verts, double *xp, double *yp,
+                  int nr_points, double *x, double *y,
+                  unsigned char *result)
+
+
+def points_inside_poly(points, verts):
+     """Test whether points lie inside a polygon.
+
+     Parameters
+     ----------
+     points : (N, 2) array
+         Input points, ``(x, y)``.
+     verts : (M, 2) array
+         Vertices of the polygon, sorted either clockwise or anti-clockwise.
+
+     Returns
+     -------
+     mask : (N,) array of bool
+         True if corresponding point is inside the polygon.
+
+     """
+     cdef np.ndarray[np.double_t, ndim=1, mode="c"] x, y, vx, vy
+
+     points = np.asarray(points)
+     verts = np.asarray(verts)
+
+     x = points[:, 0].astype(np.double)
+     y = points[:, 1].astype(np.double)
+
+     vx = verts[:, 0].astype(np.double)
+     vy = verts[:, 1].astype(np.double)
+
+     cdef np.ndarray[np.uint8_t, ndim=1] out = \
+          np.zeros(x.shape[0], dtype=np.uint8)
+     
+     npnpoly(vx.shape[0], <double*>vx.data, <double*>vy.data,
+             x.shape[0], <double*>x.data, <double*>y.data,
+             <unsigned char*>out.data)
+
+     return out.astype(bool)
+
diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py
new file mode 100644
index 00000000..36c142aa
--- /dev/null
+++ b/skimage/morphology/convex_hull.py
@@ -0,0 +1,56 @@
+__all__ = ['convex_hull']
+
+import numpy as np
+from scipy.spatial import Delaunay
+from ._pnpoly import points_inside_poly
+
+def convex_hull(image):
+    """Compute the convex hull of a binary image.
+
+    The convex hull is the set of pixels included in the smallest convex
+    polygon that surround all white pixels in the input image.
+
+    Parameters
+    ----------
+    image : ndarray
+        Binary input image.  This array is cast to bool before processing.
+
+    Returns
+    -------
+    hull : ndarray of uint8
+        Binary image with pixels in convex hull set to 255.
+
+    References
+    ----------
+    .. [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
+
+    # Add a vertex for the middle of each pixel edge
+    coords_corners = np.empty((0, 2))
+    for x_offset, y_offset in zip((0, 0, -0.5, 0.5),
+                                  (-0.5, 0.5, 0, 0)):
+            coords_corners = np.vstack((coords_corners,
+                                        coords + [x_offset, y_offset]))
+
+    coords = coords_corners
+
+    # Find the convex hull
+    chull = Delaunay(coords).convex_hull
+    v = coords[np.unique(chull)]
+
+    # Sort vertices clock-wise
+    v_centred = v - v.mean(axis=0)
+    angles = np.arctan2(v_centred[:, 0], v_centred[:, 1])
+    v = v[np.argsort(angles)]
+
+    # 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])
+    
+    return mask
diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py
index 46010989..2b3f605c 100644
--- a/skimage/morphology/setup.py
+++ b/skimage/morphology/setup.py
@@ -15,6 +15,7 @@ def configuration(parent_package='', top_path=None):
     cython(['cmorph.pyx'], working_path=base_path)
     cython(['_watershed.pyx'], working_path=base_path)
     cython(['_skeletonize.pyx'], working_path=base_path)
+    cython(['_pnpoly.pyx'], working_path=base_path)
 
     config.add_extension('ccomp', sources=['ccomp.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -24,6 +25,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_skeletonize', sources=['_skeletonize.c'],
                          include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('_pnpoly', sources=['_pnpoly.c'],
+                         include_dirs=[get_numpy_include_dirs()])
 
 
     return config
diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py
new file mode 100644
index 00000000..f6abe58f
--- /dev/null
+++ b/skimage/morphology/tests/test_convex_hull.py
@@ -0,0 +1,25 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+from skimage.morphology import convex_hull
+
+def test_basic():
+    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, 0, 0, 0, 1, 0, 0],
+         [0, 1, 0, 0, 0, 0, 0, 1, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)
+
+    expected = np.array(
+        [[0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 1, 0, 0, 0, 0],
+         [0, 0, 0, 1, 1, 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=bool)
+
+    assert_array_equal(convex_hull(image), expected)
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()
diff --git a/skimage/morphology/tests/test_pnpoly.py b/skimage/morphology/tests/test_pnpoly.py
new file mode 100644
index 00000000..eb8dbf87
--- /dev/null
+++ b/skimage/morphology/tests/test_pnpoly.py
@@ -0,0 +1,26 @@
+import numpy as np
+
+from skimage.morphology._pnpoly import points_inside_poly
+
+class test_poly():
+    def test_square(self):
+        v = np.array([[0, 0],
+                      [0, 1],
+                      [1, 1],
+                      [1, 0]])
+        assert(points_inside_poly([[0.5, 0.5]], v)[0])
+        assert(not points_inside_poly([[-0.1, 0.1]], v)[0])
+
+    def test_triangle(self):
+        v = np.array([[0, 0],
+                      [1, 0],
+                      [0.5, 0.75]])
+        assert(points_inside_poly([[0.5, 0.7]], v)[0])
+        assert(not points_inside_poly([[0.5, 0.76]], v)[0])
+        assert(not points_inside_poly([[0.7, 0.5]], v)[0])
+
+    def test_type(self):
+        assert(points_inside_poly([[0, 0]], [[0, 0]]).dtype == np.bool)
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()

From 6bbeda04f3d732e705100e2844670108db109c96 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Wed, 26 Oct 2011 19:37:34 -0700
Subject: [PATCH 2/7] ENH: Pre-allocate coordinate corner array.

---
 skimage/morphology/convex_hull.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py
index 36c142aa..d7c35f97 100644
--- a/skimage/morphology/convex_hull.py
+++ b/skimage/morphology/convex_hull.py
@@ -28,13 +28,13 @@ def convex_hull(image):
     image = image.astype(bool)
     r, c = np.nonzero(image)
     coords = np.vstack((r, c)).T
+    N = len(coords)
 
     # Add a vertex for the middle of each pixel edge
-    coords_corners = np.empty((0, 2))
-    for x_offset, y_offset in zip((0, 0, -0.5, 0.5),
-                                  (-0.5, 0.5, 0, 0)):
-            coords_corners = np.vstack((coords_corners,
-                                        coords + [x_offset, y_offset]))
+    coords_corners = np.empty((N * 4, 2))
+    for i, (x_offset, y_offset) in enumerate(zip((0, 0, -0.5, 0.5),
+                                                 (-0.5, 0.5, 0, 0))):
+        coords_corners[i * N:(i + 1) * N] = coords + [x_offset, y_offset]
 
     coords = coords_corners
 

From 74f7e01e42ad0e49a3900965a2e6d5f5046dcb03 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Wed, 26 Oct 2011 19:55:46 -0700
Subject: [PATCH 3/7] ENH: Improve convex_hull execution speed.

---
 skimage/morphology/_pnpoly.h            |  5 +--
 skimage/morphology/_pnpoly.pyx          | 46 +++++++++++++++++++++++++
 skimage/morphology/convex_hull.py       |  6 ++--
 skimage/morphology/tests/test_pnpoly.py | 16 +++++++--
 4 files changed, 65 insertions(+), 8 deletions(-)

diff --git a/skimage/morphology/_pnpoly.h b/skimage/morphology/_pnpoly.h
index a477a70f..95c89bcb 100644
--- a/skimage/morphology/_pnpoly.h
+++ b/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]))) &&
diff --git a/skimage/morphology/_pnpoly.pyx b/skimage/morphology/_pnpoly.pyx
index 65c9e55e..c492eb0d 100644
--- a/skimage/morphology/_pnpoly.pyx
+++ b/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
      -------
diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py
index d7c35f97..3101e542 100644
--- a/skimage/morphology/convex_hull.py
+++ b/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
diff --git a/skimage/morphology/tests/test_pnpoly.py b/skimage/morphology/tests/test_pnpoly.py
index eb8dbf87..33efe15f 100644
--- a/skimage/morphology/tests/test_pnpoly.py
+++ b/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()

From b34b9d8ddd09619dcb97deda44234165acdaacc9 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Wed, 26 Oct 2011 23:17:00 -0700
Subject: [PATCH 4/7] ENH: Only perform convex hull on select candidate pixels.
 Gets us a factor 30 speedup.

---
 skimage/morphology/_convex_hull.pyx          | 59 ++++++++++++++++++++
 skimage/morphology/convex_hull.py            | 14 +++--
 skimage/morphology/setup.py                  |  4 +-
 skimage/morphology/tests/test_convex_hull.py | 34 +++++++++++
 4 files changed, 106 insertions(+), 5 deletions(-)
 create mode 100644 skimage/morphology/_convex_hull.pyx

diff --git a/skimage/morphology/_convex_hull.pyx b/skimage/morphology/_convex_hull.pyx
new file mode 100644
index 00000000..dc426900
--- /dev/null
+++ b/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]
diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py
index 3101e542..034acd0b 100644
--- a/skimage/morphology/convex_hull.py
+++ b/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),
diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py
index 2b3f605c..e28446bd 100644
--- a/skimage/morphology/setup.py
+++ b/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
 
diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py
index f6abe58f..ad063458 100644
--- a/skimage/morphology/tests/test_convex_hull.py
+++ b/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()

From 2191d296113438d431c47a8d864ba7ba8b6b30d2 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Fri, 28 Oct 2011 14:36:05 -0700
Subject: [PATCH 5/7] ENH: In test, do not compute hull twice.

---
 skimage/morphology/tests/test_convex_hull.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py
index ad063458..2d150abd 100644
--- a/skimage/morphology/tests/test_convex_hull.py
+++ b/skimage/morphology/tests/test_convex_hull.py
@@ -53,7 +53,7 @@ def test_possible_hull():
                          [4, 6]])
 
     ph = possible_hull(image)
-    assert_array_equal(possible_hull(image), expected)
+    assert_array_equal(ph, expected)
 
 if __name__ == "__main__":
     np.testing.run_module_suite()

From 33469fa69eb53143ae09bd710ea46f8f05f42090 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Fri, 28 Oct 2011 14:56:41 -0700
Subject: [PATCH 6/7] BUG: Allow test suite to pass with scipy < 0.9.

---
 skimage/morphology/convex_hull.py            |  7 ++++++-
 skimage/morphology/tests/test_convex_hull.py | 10 +++++++++-
 2 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py
index 034acd0b..487a547b 100644
--- a/skimage/morphology/convex_hull.py
+++ b/skimage/morphology/convex_hull.py
@@ -1,7 +1,6 @@
 __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
 
@@ -44,6 +43,12 @@ def convex_hull(image):
 
     coords = coords_corners
 
+    try:
+        from scipy.spatial import Delaunay
+    except ImportError:
+        raise ImportError('Could not import scipy.spatial, only available in '
+                          'scipy >= 0.9.')
+
     # Find the convex hull
     chull = Delaunay(coords).convex_hull
     v = coords[np.unique(chull)]
diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py
index 2d150abd..e8fb5514 100644
--- a/skimage/morphology/tests/test_convex_hull.py
+++ b/skimage/morphology/tests/test_convex_hull.py
@@ -1,8 +1,16 @@
 import numpy as np
 from numpy.testing import assert_array_equal
+from numpy.testing.decorators import skipif
 from skimage.morphology import convex_hull
 from skimage.morphology._convex_hull import possible_hull
 
+try:
+    import scipy.spatial
+    scipy_spatial = True
+except ImportError:
+    scipy_spatial = False
+
+@skipif(not scipy_spatial)
 def test_basic():
     image = np.array(
         [[0, 0, 0, 0, 0, 0, 0, 0, 0],
@@ -22,7 +30,7 @@ def test_basic():
 
     assert_array_equal(convex_hull(image), expected)
 
-
+@skipif(not scipy_spatial)
 def test_possible_hull():
     image = np.array(
         [[0, 0, 0, 0, 0, 0, 0, 0, 0],

From 349576735882e1774fbbef9d9a8f97af8213d97f Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Wed, 2 Nov 2011 09:38:54 -0700
Subject: [PATCH 7/7] ENH: Rename convex_hull to convex_hull_image to allow
 later clean API addition of convex_hull with labels.

---
 skimage/morphology/__init__.py               | 2 +-
 skimage/morphology/convex_hull.py            | 6 +++---
 skimage/morphology/tests/test_convex_hull.py | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py
index 5a7d67d7..e3176579 100644
--- a/skimage/morphology/__init__.py
+++ b/skimage/morphology/__init__.py
@@ -3,4 +3,4 @@ from selem import *
 from .ccomp import label
 from watershed import watershed, is_local_maximum
 from skeletonize import skeletonize, medial_axis
-from .convex_hull import convex_hull
+from .convex_hull import convex_hull_image
diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py
index 487a547b..f461dc4d 100644
--- a/skimage/morphology/convex_hull.py
+++ b/skimage/morphology/convex_hull.py
@@ -1,11 +1,11 @@
-__all__ = ['convex_hull']
+__all__ = ['convex_hull_image']
 
 import numpy as np
 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.
+def convex_hull_image(image):
+    """Compute the convex hull image of a binary image.
 
     The convex hull is the set of pixels included in the smallest convex
     polygon that surround all white pixels in the input image.
diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py
index e8fb5514..21f45cd7 100644
--- a/skimage/morphology/tests/test_convex_hull.py
+++ b/skimage/morphology/tests/test_convex_hull.py
@@ -1,7 +1,7 @@
 import numpy as np
 from numpy.testing import assert_array_equal
 from numpy.testing.decorators import skipif
-from skimage.morphology import convex_hull
+from skimage.morphology import convex_hull_image
 from skimage.morphology._convex_hull import possible_hull
 
 try:
@@ -28,7 +28,7 @@ def test_basic():
          [0, 1, 1, 1, 1, 1, 1, 1, 0],
          [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)
 
-    assert_array_equal(convex_hull(image), expected)
+    assert_array_equal(convex_hull_image(image), expected)
 
 @skipif(not scipy_spatial)
 def test_possible_hull():