ENH: Add binary convex hull computation.

CONTRIBUTORS.txt

@@ -84,3 +84,5 @@
 - Nelle Varoquaux
   Renaming of the package to ``skimage``.
 
+- W. Randolph Franklin
+  Point in polygon test.

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()

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

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

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)
+

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

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

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()

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()