From b34b9d8ddd09619dcb97deda44234165acdaacc9 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Wed, 26 Oct 2011 23:17:00 -0700 Subject: [PATCH] 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()