From 567b0fc54f66130d09ce0ed4181e3865ba6c6b13 Mon Sep 17 00:00:00 2001 From: Egor Panfilov Date: Sun, 20 Mar 2016 23:37:43 +0300 Subject: [PATCH 1/2] Created Python wrapper, renamed files, fixed error message --- bento.info | 4 +- skimage/segmentation/_quickshift.py | 74 +++++++++++++++++++ .../{_quickshift.pyx => _quickshift_cy.pyx} | 71 +++++------------- skimage/segmentation/setup.py | 4 +- skimage/segmentation/tests/test_quickshift.py | 2 +- 5 files changed, 98 insertions(+), 57 deletions(-) create mode 100644 skimage/segmentation/_quickshift.py rename skimage/segmentation/{_quickshift.pyx => _quickshift_cy.pyx} (66%) diff --git a/bento.info b/bento.info index 82ffceb2..0dad510e 100644 --- a/bento.info +++ b/bento.info @@ -113,9 +113,9 @@ Library: Extension: skimage.segmentation._slic Sources: skimage/segmentation/_slic.pyx - Extension: skimage.segmentation._quickshift + Extension: skimage.segmentation._quickshift_cy Sources: - skimage/segmentation/_quickshift.pyx + skimage/segmentation/_quickshift_cy.pyx Extension: skimage.morphology._skeletonize_cy Sources: skimage/morphology/_skeletonize_cy.pyx diff --git a/skimage/segmentation/_quickshift.py b/skimage/segmentation/_quickshift.py new file mode 100644 index 00000000..d5b11e78 --- /dev/null +++ b/skimage/segmentation/_quickshift.py @@ -0,0 +1,74 @@ +import numpy as np + +import scipy.ndimage as ndi + +from ..util import img_as_float +from ..color import rgb2lab + +from ._quickshift_cy import _quickshift_cython + + +def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10, + return_tree=False, sigma=0, convert2lab=True, random_seed=42): + """Segments image using quickshift clustering in Color-(x,y) space. + + Produces an oversegmentation of the image using the quickshift mode-seeking + algorithm. + + Parameters + ---------- + image : (width, height, channels) ndarray + Input image. + ratio : float, optional, between 0 and 1 (default 1). + Balances color-space proximity and image-space proximity. + Higher values give more weight to color-space. + kernel_size : float, optional (default 5) + Width of Gaussian kernel used in smoothing the + sample density. Higher means fewer clusters. + max_dist : float, optional (default 10) + Cut-off point for data distances. + Higher means fewer clusters. + return_tree : bool, optional (default False) + Whether to return the full segmentation hierarchy tree and distances. + sigma : float, optional (default 0) + Width for Gaussian smoothing as preprocessing. Zero means no smoothing. + convert2lab : bool, optional (default True) + Whether the input should be converted to Lab colorspace prior to + segmentation. For this purpose, the input is assumed to be RGB. + random_seed : int, optional (default 42) + Random seed used for breaking ties. + + Returns + ------- + segment_mask : (width, height) ndarray + Integer mask indicating segment labels. + + Notes + ----- + The authors advocate to convert the image to Lab color space prior to + segmentation, though this is not strictly necessary. For this to work, the + image must be given in RGB format. + + References + ---------- + .. [1] Quick shift and kernel methods for mode seeking, + Vedaldi, A. and Soatto, S. + European Conference on Computer Vision, 2008 + """ + + image = img_as_float(np.atleast_3d(image)) + if convert2lab: + if image.shape[2] != 3: + ValueError("Only RGB images can be converted to Lab space.") + image = rgb2lab(image) + + if kernel_size < 1: + raise ValueError("`kernel_size` should be >= 1.") + + image = ndi.gaussian_filter(image, [sigma, sigma, 0]) + image = np.ascontiguousarray(image * ratio) + + segment_mask = _quickshift_cython( + image, kernel_size=kernel_size, max_dist=max_dist, + return_tree=return_tree, random_seed=random_seed) + return segment_mask diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift_cy.pyx similarity index 66% rename from skimage/segmentation/_quickshift.pyx rename to skimage/segmentation/_quickshift_cy.pyx index 3ec70d98..5a78ccc7 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift_cy.pyx @@ -2,20 +2,19 @@ #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False -import numpy as np -from scipy import ndimage as ndi -from itertools import product +import numpy as np cimport numpy as cnp + from libc.math cimport exp, sqrt from libc.float cimport DBL_MAX -from ..util import img_as_float -from ..color import rgb2lab - -def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10, - return_tree=False, sigma=0, convert2lab=True, random_seed=None): +def _quickshift_cython(cnp.ndarray[double, ndim=3, mode="c"] image, + float kernel_size, + float max_dist, + bint return_tree, + int random_seed): """Segments image using quickshift clustering in Color-(x,y) space. Produces an oversegmentation of the image using the quickshift mode-seeking @@ -25,53 +24,22 @@ def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10, ---------- image : (width, height, channels) ndarray Input image. - ratio : float, optional, between 0 and 1 (default 1). - Balances color-space proximity and image-space proximity. - Higher values give more weight to color-space. - kernel_size : float, optional (default 5) + kernel_size : float Width of Gaussian kernel used in smoothing the sample density. Higher means fewer clusters. - max_dist : float, optional (default 10) + max_dist : float Cut-off point for data distances. Higher means fewer clusters. - return_tree : bool, optional (default False) + return_tree : bool Whether to return the full segmentation hierarchy tree and distances. - sigma : float, optional (default 0) - Width for Gaussian smoothing as preprocessing. Zero means no smoothing. - convert2lab : bool, optional (default True) - Whether the input should be converted to Lab colorspace prior to - segmentation. For this purpose, the input is assumed to be RGB. - random_seed : None (default) or int, optional + random_seed : int Random seed used for breaking ties. Returns ------- segment_mask : (width, height) ndarray Integer mask indicating segment labels. - - Notes - ----- - The authors advocate to convert the image to Lab color space prior to - segmentation, though this is not strictly necessary. For this to work, the - image must be given in RGB format. - - References - ---------- - .. [1] Quick shift and kernel methods for mode seeking, - Vedaldi, A. and Soatto, S. - European Conference on Computer Vision, 2008 - - """ - image = img_as_float(np.atleast_3d(image)) - if convert2lab: - if image.shape[2] != 3: - ValueError("Only RGB images can be converted to Lab space.") - image = rgb2lab(image) - - image = ndi.gaussian_filter(img_as_float(image), [sigma, sigma, 0]) - cdef cnp.ndarray[dtype=cnp.float_t, ndim=3, mode="c"] image_c \ - = np.ascontiguousarray(image) * ratio random_state = np.random.RandomState(random_seed) @@ -82,26 +50,25 @@ def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10, # an effect for very high max_dist. # window size for neighboring pixels to consider - if kernel_size < 1: - raise ValueError("Sigma should be >= 1") cdef float kernel_size_sq = kernel_size**2 cdef int w = np.ceil(3 * kernel_size) - cdef Py_ssize_t height = image_c.shape[0] - cdef Py_ssize_t width = image_c.shape[1] - cdef Py_ssize_t channels = image_c.shape[2] + cdef Py_ssize_t height = image.shape[0] + cdef Py_ssize_t width = image.shape[1] + cdef Py_ssize_t channels = image.shape[2] cdef double current_density, closest, dist cdef Py_ssize_t r, c, r_, c_, channel, r_min, c_min - cdef cnp.float_t* image_p = image_c.data - cdef cnp.float_t* current_pixel_p = image_p + cdef cnp.float_t* image_p = image.data + cdef cnp.float_t* current_pixel_p cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] densities \ = np.zeros((height, width)) # compute densities with nogil: + current_pixel_p = image_p for r in range(height): for c in range(width): r_min, r_max = max(r - w, 0), min(r + w + 1, height) @@ -111,7 +78,7 @@ def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10, dist = 0 for channel in range(channels): dist += (current_pixel_p[channel] - - image_c[r_, c_, channel])**2 + image[r_, c_, channel])**2 dist += (r - r_)**2 + (c - c_)**2 densities[r, c] += exp(-dist / (2 * kernel_size_sq)) current_pixel_p += channels @@ -143,7 +110,7 @@ def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10, # (width * height * windowsize**2) for channel in range(channels): dist += (current_pixel_p[channel] - - image_c[r_, c_, channel])**2 + image[r_, c_, channel])**2 dist += (r - r_)**2 + (c - c_)**2 if dist < closest: closest = dist diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index ec4ac4ff..3f01a146 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -14,8 +14,8 @@ def configuration(parent_package='', top_path=None): cython(['_felzenszwalb_cy.pyx'], working_path=base_path) config.add_extension('_felzenszwalb_cy', sources=['_felzenszwalb_cy.c'], include_dirs=[get_numpy_include_dirs()]) - cython(['_quickshift.pyx'], working_path=base_path) - config.add_extension('_quickshift', sources=['_quickshift.c'], + cython(['_quickshift_cy.pyx'], working_path=base_path) + config.add_extension('_quickshift_cy', sources=['_quickshift_cy.c'], include_dirs=[get_numpy_include_dirs()]) cython(['_slic.pyx'], working_path=base_path) config.add_extension('_slic', sources=['_slic.c'], diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py index 21dbe8ea..ba21532c 100644 --- a/skimage/segmentation/tests/test_quickshift.py +++ b/skimage/segmentation/tests/test_quickshift.py @@ -41,7 +41,7 @@ def test_color(): assert_array_equal(seg[10:, 10:], 3) seg2 = quickshift(img, kernel_size=1, max_dist=2, random_seed=0, - convert2lab=False, sigma=0) + convert2lab=False, sigma=0) # very oversegmented: assert_equal(len(np.unique(seg2)), 7) # still don't cross lines From d19817cacd53a6a272ec3e369123fe6d1b14924a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 22 Mar 2016 21:41:28 +0100 Subject: [PATCH 2/2] Use typed memoryviews for quickshift and other Cython improvements --- skimage/segmentation/_quickshift_cy.pyx | 85 +++++++++++++------------ 1 file changed, 43 insertions(+), 42 deletions(-) diff --git a/skimage/segmentation/_quickshift_cy.pyx b/skimage/segmentation/_quickshift_cy.pyx index 5a78ccc7..7b810303 100644 --- a/skimage/segmentation/_quickshift_cy.pyx +++ b/skimage/segmentation/_quickshift_cy.pyx @@ -6,15 +6,12 @@ import numpy as np cimport numpy as cnp -from libc.math cimport exp, sqrt +from libc.math cimport exp, sqrt, ceil from libc.float cimport DBL_MAX -def _quickshift_cython(cnp.ndarray[double, ndim=3, mode="c"] image, - float kernel_size, - float max_dist, - bint return_tree, - int random_seed): +def _quickshift_cython(double[:, :, ::1] image, double kernel_size, + double max_dist, bint return_tree, int random_seed): """Segments image using quickshift clustering in Color-(x,y) space. Produces an oversegmentation of the image using the quickshift mode-seeking @@ -50,57 +47,57 @@ def _quickshift_cython(cnp.ndarray[double, ndim=3, mode="c"] image, # an effect for very high max_dist. # window size for neighboring pixels to consider - cdef float kernel_size_sq = kernel_size**2 - cdef int w = np.ceil(3 * kernel_size) + cdef double inv_kernel_size_sqr = -0.5 / kernel_size**2 + cdef int kernel_width = ceil(3 * kernel_size) cdef Py_ssize_t height = image.shape[0] cdef Py_ssize_t width = image.shape[1] cdef Py_ssize_t channels = image.shape[2] + + cdef double[:, ::1] densities = np.zeros((height, width), dtype=np.double) + cdef double current_density, closest, dist - - cdef Py_ssize_t r, c, r_, c_, channel, r_min, c_min - - cdef cnp.float_t* image_p = image.data - cdef cnp.float_t* current_pixel_p - - cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] densities \ - = np.zeros((height, width)) + cdef Py_ssize_t r, c, r_, c_, channel, r_min, r_max, c_min, c_max + cdef double* current_pixel_ptr # compute densities with nogil: - current_pixel_p = image_p + current_pixel_ptr = &image[0, 0, 0] for r in range(height): + r_min = max(r - kernel_width, 0) + r_max = min(r + kernel_width + 1, height) for c in range(width): - r_min, r_max = max(r - w, 0), min(r + w + 1, height) - c_min, c_max = max(c - w, 0), min(c + w + 1, width) + c_min = max(c - kernel_width, 0) + c_max = min(c + kernel_width + 1, width) for r_ in range(r_min, r_max): for c_ in range(c_min, c_max): dist = 0 for channel in range(channels): - dist += (current_pixel_p[channel] - + dist += (current_pixel_ptr[channel] - image[r_, c_, channel])**2 dist += (r - r_)**2 + (c - c_)**2 - densities[r, c] += exp(-dist / (2 * kernel_size_sq)) - current_pixel_p += channels + densities[r, c] += exp(dist * inv_kernel_size_sqr) + current_pixel_ptr += channels # this will break ties that otherwise would give us headache densities += random_state.normal(scale=0.00001, size=(height, width)) # default parent to self - cdef cnp.ndarray[dtype=cnp.int_t, ndim=2] parent \ - = np.arange(width * height).reshape(height, width) - cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] dist_parent \ - = np.zeros((height, width)) + cdef Py_ssize_t[:, ::1] parent = \ + np.arange(width * height, dtype=np.intp).reshape(height, width) + cdef double[:, ::1] dist_parent = np.zeros((height, width), dtype=np.double) # find nearest node with higher density with nogil: - current_pixel_p = image_p + current_pixel_ptr = &image[0, 0, 0] for r in range(height): + r_min = max(r - kernel_width, 0) + r_max = min(r + kernel_width + 1, height) for c in range(width): current_density = densities[r, c] closest = DBL_MAX - r_min, r_max = max(r - w, 0), min(r + w + 1, height) - c_min, c_max = max(c - w, 0), min(c + w + 1, width) + c_min = max(c - kernel_width, 0) + c_max = min(c + kernel_width + 1, width) for r_ in range(r_min, r_max): for c_ in range(c_min, c_max): if densities[r_, c_] > current_density: @@ -109,27 +106,31 @@ def _quickshift_cython(cnp.ndarray[double, ndim=3, mode="c"] image, # we get crazy memory overhead # (width * height * windowsize**2) for channel in range(channels): - dist += (current_pixel_p[channel] - + dist += (current_pixel_ptr[channel] - image[r_, c_, channel])**2 dist += (r - r_)**2 + (c - c_)**2 if dist < closest: closest = dist parent[r, c] = r_ * width + c_ dist_parent[r, c] = sqrt(closest) - current_pixel_p += channels + current_pixel_ptr += channels + + dist_parent_flat = np.array(dist_parent).ravel() + parent_flat = np.array(parent).ravel() - dist_parent_flat = dist_parent.ravel() - flat = parent.ravel() # remove parents with distance > max_dist too_far = dist_parent_flat > max_dist - flat[too_far] = np.arange(width * height)[too_far] - old = np.zeros_like(flat) + parent_flat[too_far] = np.arange(width * height)[too_far] + old = np.zeros_like(parent_flat) + # flatten forest (mark each pixel with root of corresponding tree) - while (old != flat).any(): - old = flat - flat = flat[flat] - flat = np.unique(flat, return_inverse=True)[1] - flat = flat.reshape(height, width) + while (old != parent_flat).any(): + old = parent_flat + parent_flat = parent_flat[parent_flat] + + parent_flat = np.unique(parent_flat, return_inverse=True)[1] + parent_flat = parent_flat.reshape(height, width) + if return_tree: - return flat, parent, dist_parent - return flat + return parent_flat, parent, dist_parent + return parent_flat