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.pyx deleted file mode 100644 index 3ec70d98..00000000 --- a/skimage/segmentation/_quickshift.pyx +++ /dev/null @@ -1,168 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False -import numpy as np -from scipy import ndimage as ndi -from itertools import product - -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): - """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 : None (default) or int, optional - 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) - - # TODO join orphaned roots? - # Some nodes might not have a point of higher density within the - # search window. We could do a global search over these in the end. - # Reference implementation doesn't do that, though, and it only has - # 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 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.ndarray[dtype=cnp.float_t, ndim=2] densities \ - = np.zeros((height, width)) - - # compute densities - with nogil: - for r in range(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) - 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] - - image_c[r_, c_, channel])**2 - dist += (r - r_)**2 + (c - c_)**2 - densities[r, c] += exp(-dist / (2 * kernel_size_sq)) - current_pixel_p += 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)) - - # find nearest node with higher density - with nogil: - current_pixel_p = image_p - for r in range(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) - for r_ in range(r_min, r_max): - for c_ in range(c_min, c_max): - if densities[r_, c_] > current_density: - dist = 0 - # We compute the distances twice since otherwise - # we get crazy memory overhead - # (width * height * windowsize**2) - for channel in range(channels): - dist += (current_pixel_p[channel] - - image_c[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 - - 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) - # 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) - if return_tree: - return flat, parent, dist_parent - return flat diff --git a/skimage/segmentation/_quickshift_cy.pyx b/skimage/segmentation/_quickshift_cy.pyx new file mode 100644 index 00000000..7b810303 --- /dev/null +++ b/skimage/segmentation/_quickshift_cy.pyx @@ -0,0 +1,136 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as cnp + +from libc.math cimport exp, sqrt, ceil +from libc.float cimport DBL_MAX + + +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 + algorithm. + + Parameters + ---------- + image : (width, height, channels) ndarray + Input image. + kernel_size : float + Width of Gaussian kernel used in smoothing the + sample density. Higher means fewer clusters. + max_dist : float + Cut-off point for data distances. + Higher means fewer clusters. + return_tree : bool + Whether to return the full segmentation hierarchy tree and distances. + random_seed : int + Random seed used for breaking ties. + + Returns + ------- + segment_mask : (width, height) ndarray + Integer mask indicating segment labels. + """ + + random_state = np.random.RandomState(random_seed) + + # TODO join orphaned roots? + # Some nodes might not have a point of higher density within the + # search window. We could do a global search over these in the end. + # Reference implementation doesn't do that, though, and it only has + # an effect for very high max_dist. + + # window size for neighboring pixels to consider + 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, r_max, c_min, c_max + cdef double* current_pixel_ptr + + # compute densities + with nogil: + 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): + 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_ptr[channel] - + image[r_, c_, channel])**2 + dist += (r - r_)**2 + (c - c_)**2 + 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 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_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 + 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: + dist = 0 + # We compute the distances twice since otherwise + # we get crazy memory overhead + # (width * height * windowsize**2) + for channel in range(channels): + 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_ptr += channels + + dist_parent_flat = np.array(dist_parent).ravel() + parent_flat = np.array(parent).ravel() + + # remove parents with distance > max_dist + too_far = dist_parent_flat > max_dist + 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 != 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 parent_flat, parent, dist_parent + return parent_flat 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