Merge pull request #2021 from ahojnnes/soupault-fix_2005

skimage.segmentation.quickshift signature is missing from API docs, third attempt
This commit is contained in:
Steven Silvester
2016-03-23 06:22:13 -05:00
6 changed files with 215 additions and 173 deletions
+2 -2
View File
@@ -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
+74
View File
@@ -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
-168
View File
@@ -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 = <cnp.float_t*> 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
+136
View File
@@ -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 = <int>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
+2 -2
View File
@@ -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'],
@@ -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