Merge pull request #672 from jni/chull-fix

Convex hull fix

skimage/morphology/convex_hull.py

@@ -4,6 +4,7 @@ import numpy as np
 from ._pnpoly import grid_points_inside_poly
 from ._convex_hull import possible_hull
 from skimage.morphology import label
+from skimage.util import unique_rows
 
 
 def convex_hull_image(image):
@@ -43,7 +44,9 @@ def convex_hull_image(image):
                                                  (-0.5, 0.5, 0, 0))):
         coords_corners[i * N:(i + 1) * N] = coords + [x_offset, y_offset]
 
-    coords = coords_corners
+    # repeated coordinates can *sometimes* cause problems in 
+    # scipy.spatial.Delaunay, so we remove them.
+    coords = unique_rows(coords_corners)
 
     try:
         from scipy.spatial import Delaunay

skimage/morphology/tests/test_convex_hull.py

@@ -32,6 +32,19 @@ def test_basic():
     assert_array_equal(convex_hull_image(image), expected)
 
 
+@skipif(not scipy_spatial)
+def test_pathological_qhull_example():
+    image = np.array(
+                [[0, 0, 0, 0, 1, 0, 0],
+                 [0, 0, 1, 1, 1, 1, 1],
+                 [1, 1, 1, 0, 0, 0, 0]], dtype=bool)
+    expected = np.array(
+                [[0, 0, 0, 1, 1, 1, 0],
+                 [0, 1, 1, 1, 1, 1, 1],
+                 [1, 1, 1, 1, 0, 0, 0]], dtype=bool)
+    assert_array_equal(convex_hull_image(image), expected)
+
+
 @skipif(not scipy_spatial)
 def test_possible_hull():
     image = np.array(

skimage/util/__init__.py

@@ -12,6 +12,7 @@ else:
     from numpy import pad
 del numpy, ver, chk
 from ._regular_grid import regular_grid
+from .unique import unique_rows
 
 
 __all__ = ['img_as_float',
@@ -24,4 +25,5 @@ __all__ = ['img_as_float',
            'view_as_windows',
            'pad',
            'random_noise',
-           'regular_grid']
+           'regular_grid',
+           'unique_rows']

skimage/util/tests/test_unique_rows.py

@@ -0,0 +1,40 @@
+import numpy as np
+from numpy.testing import assert_equal, assert_raises
+from skimage.util import unique_rows
+
+
+def test_discontiguous_array():
+    ar = np.array([[1, 0, 1], [0, 1, 0], [1, 0, 1]], np.uint8)
+    ar = ar[::2]
+    ar_out = unique_rows(ar)
+    desired_ar_out = np.array([[1, 0, 1]], np.uint8)
+    assert_equal(ar_out, desired_ar_out)
+
+
+def test_uint8_array():
+    ar = np.array([[1, 0, 1], [0, 1, 0], [1, 0, 1]], np.uint8)
+    ar_out = unique_rows(ar)
+    desired_ar_out = np.array([[0, 1, 0], [1, 0, 1]], np.uint8)
+    assert_equal(ar_out, desired_ar_out)
+
+
+def test_float_array():
+    ar = np.array([[1.1, 0.0, 1.1], [0.0, 1.1, 0.0], [1.1, 0.0, 1.1]],
+                  np.float)
+    ar_out = unique_rows(ar)
+    desired_ar_out = np.array([[0.0, 1.1, 0.0], [1.1, 0.0, 1.1]], np.float)
+    assert_equal(ar_out, desired_ar_out)
+
+
+def test_1d_array():
+    ar = np.array([1, 0, 1, 1], np.uint8)
+    assert_raises(ValueError, unique_rows, ar)
+
+
+def test_3d_array():
+    ar = np.arange(8).reshape((2, 2, 2))
+    assert_raises(ValueError, unique_rows, ar)
+
+
+if __name__ == '__main__':
+    np.testing.run_module_suite()

skimage/util/unique.py

@@ -0,0 +1,50 @@
+import numpy as np
+
+
+def unique_rows(ar):
+    """Remove repeated rows from a 2D array.
+
+    In particular, if given an array of coordinates of shape
+    (Npoints, Ndim), it will remove repeated points.
+
+    Parameters
+    ----------
+    ar : 2D np.ndarray
+        The input array.
+
+    Returns
+    -------
+    ar_out : 2D np.ndarray
+        A copy of the input array with repeated rows removed.
+
+    Raises
+    ------
+    ValueError : if `ar` is not two-dimensional.
+
+    Notes
+    -----
+    The function will generate a copy of `ar` if it is not
+    C-contiguous, which will negatively affect performance for large
+    input arrays.
+
+    Examples
+    --------
+    >>> ar = np.array([[1, 0, 1],
+                       [0, 1, 0],
+                       [1, 0, 1]], np.uint8)
+    >>> aru = unique_rows(ar)
+    array([[0, 1, 0],
+           [1, 0, 1]], dtype=uint8)
+    """
+    if ar.ndim != 2:
+        raise ValueError("unique_rows() only makes sense for 2D arrays, "
+                         "got %dd" % ar.ndim)
+    # the view in the next line only works if the array is C-contiguous
+    ar = np.ascontiguousarray(ar)
+    # np.unique() finds identical items in a raveled array. To make it
+    # see each row as a single item, we create a view of each row as a
+    # byte string of length itemsize times number of columns in `ar`
+    ar_row_view = ar.view('|S%d' % (ar.itemsize * ar.shape[1]))
+    _, unique_row_indices = np.unique(ar_row_view, return_index=True)
+    ar_out = ar[unique_row_indices]
+    return ar_out