mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-19 11:27:45 +08:00
Merge pull request #206 from amueller/felsenzwalb
ENH: MRG Segmentation algorithms.
This commit is contained in:
Stefan van der Walt
@@ -74,6 +74,7 @@
|
||||
|
||||
- Andreas Mueller
|
||||
Example data set loader.
|
||||
Quickshift image segmentation, Felzenszwalbs fast graph based segmentation.
|
||||
|
||||
- Yaroslav Halchenko
|
||||
For sharing his expert advice on Debian packaging.
|
||||
|
||||
@@ -0,0 +1,91 @@
|
||||
"""
|
||||
====================================================
|
||||
Comparison of segmentation and superpixel algorithms
|
||||
====================================================
|
||||
|
||||
This example compares three popular low-level image segmentation methods. As
|
||||
it is difficult to obtain good segmentations, and the definition of "good"
|
||||
often depends on the application, these methods are usually used for obtaining
|
||||
an oversegmentation, also known as superpixels. These superpixels then serve as
|
||||
a basis for more sophisticated algorithms such as CRFs.
|
||||
|
||||
|
||||
Felzenszwalb's efficient graph based segmentation
|
||||
-------------------------------------------------
|
||||
This fast 2D image segmentation algorithm, proposed in [1]_ is popular in the
|
||||
computer vision community.
|
||||
The algorithm has a single ``scale`` parameter that influences the segment
|
||||
size. The actual size and number of segments can vary greatly, depending on
|
||||
local contrast.
|
||||
|
||||
.. [1] Efficient graph-based image segmentation, Felzenszwalb, P.F. and
|
||||
Huttenlocher, D.P. International Journal of Computer Vision, 2004
|
||||
|
||||
|
||||
Quickshift image segmentation
|
||||
-----------------------------
|
||||
|
||||
Quickshift is a relatively recent 2D image segmentation algorithm, based on an
|
||||
approximation of kernelized mean-shift. Therefore it belongs to the family of
|
||||
local mode-seeking algorithms and is applied to the 5D space consisting of
|
||||
color information and image location [2]_.
|
||||
|
||||
One of the benefits of quickshift is that it actually computes a
|
||||
hierarchical segmentation on multiple scales simultaneously.
|
||||
|
||||
Quickshift has two main parameters: ``sigma`` controls the scale of the local
|
||||
density approximation, ``max_dist`` selects a level in the hierarchical
|
||||
segmentation that is produced. There is also a trade-off between distance in
|
||||
color-space and distance in image-space, given by ``ratio``.
|
||||
|
||||
.. [2] Quick shift and kernel methods for mode seeking,
|
||||
Vedaldi, A. and Soatto, S.
|
||||
European Conference on Computer Vision, 2008
|
||||
|
||||
|
||||
SLIC - K-Means based image segmentation
|
||||
---------------------------------------
|
||||
This algorithm simply performs K-means in the 5d space of color information
|
||||
and image location and is therefore closely related to quickshift. As the
|
||||
clustering method is simpler, it is very efficient. It is essential for this
|
||||
algorithm to work in Lab color space to obtain good results. The algorithm
|
||||
quickly gained momentum and is now widely used. See [3] for details. The
|
||||
``ratio`` parameter trades off color-similarity and proximity, as in the case
|
||||
of Quickshift, while ``n_segments`` chooses the number of centers for kmeans.
|
||||
|
||||
.. [3] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi,
|
||||
Pascal Fua, and Sabine Suesstrunk, SLIC Superpixels Compared to
|
||||
State-of-the-art Superpixel Methods, TPAMI, May 2012.
|
||||
"""
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
from skimage.data import lena
|
||||
from skimage.segmentation import felzenszwalb, \
|
||||
visualize_boundaries, slic, quickshift
|
||||
from skimage.util import img_as_float
|
||||
|
||||
img = img_as_float(lena()[::2, ::2])
|
||||
segments_fz = felzenszwalb(img, scale=100, sigma=0.5, min_size=50)
|
||||
segments_slic = slic(img, ratio=10, n_segments=250, sigma=1)
|
||||
segments_quick = quickshift(img, kernel_size=3, max_dist=6, ratio=0.5)
|
||||
|
||||
print("Felzenszwalb's number of segments: %d" % len(np.unique(segments_fz)))
|
||||
print("Slic number of segments: %d" % len(np.unique(segments_slic)))
|
||||
print("Quickshift number of segments: %d" % len(np.unique(segments_quick)))
|
||||
|
||||
fig, ax = plt.subplots(1, 3)
|
||||
fig.set_size_inches(8, 3, forward=True)
|
||||
plt.subplots_adjust(0.05, 0.05, 0.95, 0.95, 0.05, 0.05)
|
||||
|
||||
ax[0].imshow(visualize_boundaries(img, segments_fz))
|
||||
ax[0].set_title("Felzenszwalbs's method")
|
||||
ax[1].imshow(visualize_boundaries(img, segments_slic))
|
||||
ax[1].set_title("SLIC")
|
||||
ax[2].imshow(visualize_boundaries(img, segments_quick))
|
||||
ax[2].set_title("Quickshift")
|
||||
for a in ax:
|
||||
a.set_xticks(())
|
||||
a.set_yticks(())
|
||||
plt.show()
|
||||
@@ -0,0 +1,10 @@
|
||||
"""Export fast union find in Cython"""
|
||||
cimport numpy as np
|
||||
|
||||
DTYPE = np.int
|
||||
ctypedef np.int_t DTYPE_t
|
||||
|
||||
cdef DTYPE_t find_root(np.int_t *forest, np.int_t n)
|
||||
cdef set_root(np.int_t *forest, np.int_t n, np.int_t root)
|
||||
cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m)
|
||||
cdef link_bg(np.int_t *forest, np.int_t n, np.int_t *background_node)
|
||||
@@ -24,7 +24,6 @@ See also:
|
||||
# The term "forest" is used to indicate an array that stores one or more trees
|
||||
|
||||
DTYPE = np.int
|
||||
ctypedef np.int_t DTYPE_t
|
||||
|
||||
cdef DTYPE_t find_root(np.int_t *forest, np.int_t n):
|
||||
"""Find the root of node n.
|
||||
|
||||
@@ -1 +1,5 @@
|
||||
from .random_walker_segmentation import random_walker
|
||||
from ._felzenszwalb import felzenszwalb
|
||||
from ._slic import slic
|
||||
from ._quickshift import quickshift
|
||||
from .boundaries import find_boundaries, visualize_boundaries
|
||||
|
||||
@@ -0,0 +1,77 @@
|
||||
import warnings
|
||||
import numpy as np
|
||||
|
||||
from ._felzenszwalb_cy import _felzenszwalb_grey
|
||||
|
||||
|
||||
def felzenszwalb(image, scale=1, sigma=0.8, min_size=20):
|
||||
"""Computes Felsenszwalb's efficient graph based image segmentation.
|
||||
|
||||
Produces an oversegmentation of a multichannel (i.e. RGB) image
|
||||
using a fast, minimum spanning tree based clustering on the image grid.
|
||||
The parameter ``scale`` sets an observation level. Higher scale means
|
||||
less and larger segments. ``sigma`` is the diameter of a Gaussian kernel,
|
||||
used for smoothing the image prior to segmentation.
|
||||
|
||||
The number of produced segments as well as their size can only be
|
||||
controlled indirectly through ``scale``. Segment size within an image can
|
||||
vary greatly depending on local contrast.
|
||||
|
||||
For RGB images, the algorithm computes a separate segmentation for each
|
||||
channel and then combines these. The combined segmentation is the
|
||||
intersection of the separate segmentations on the color channels.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : (width, height, 3) or (width, height) ndarray
|
||||
Input image.
|
||||
scale : float
|
||||
Free parameter. Higher means larger clusters.
|
||||
sigma : float
|
||||
Width of Gaussian kernel used in preprocessing.
|
||||
min_size : int
|
||||
Minimum component size. Enforced using postprocessing.
|
||||
|
||||
Returns
|
||||
-------
|
||||
segment_mask : (width, height) ndarray
|
||||
Integer mask indicating segment labels.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Efficient graph-based image segmentation, Felzenszwalb, P.F. and
|
||||
Huttenlocher, D.P. International Journal of Computer Vision, 2004
|
||||
"""
|
||||
|
||||
#image = img_as_float(image)
|
||||
if image.ndim == 2:
|
||||
# assume single channel image
|
||||
return _felzenszwalb_grey(image, scale=scale, sigma=sigma)
|
||||
|
||||
elif image.ndim != 3:
|
||||
raise ValueError("Felzenswalb segmentation can only operate on RGB and"
|
||||
" grey images, but input array of ndim %d given."
|
||||
% image.ndim)
|
||||
|
||||
# assume we got 2d image with multiple channels
|
||||
n_channels = image.shape[2]
|
||||
if n_channels != 3:
|
||||
warnings.warn("Got image with %d channels. Is that really what you"
|
||||
" wanted?" % image.shape[2])
|
||||
segmentations = []
|
||||
# compute quickshift for each channel
|
||||
for c in range(n_channels):
|
||||
channel = np.ascontiguousarray(image[:, :, c])
|
||||
s = _felzenszwalb_grey(channel, scale=scale, sigma=sigma,
|
||||
min_size=min_size)
|
||||
segmentations.append(s)
|
||||
|
||||
# put pixels in same segment only if in the same segment in all images
|
||||
# we do this by combining the channels to one number
|
||||
n0 = segmentations[0].max() + 1
|
||||
n1 = segmentations[1].max() + 1
|
||||
segmentation = (segmentations[0] + segmentations[1] * n0
|
||||
+ segmentations[2] * n0 * n1)
|
||||
# make segment labels consecutive numbers starting at 0
|
||||
labels = np.unique(segmentation, return_inverse=True)[1]
|
||||
return labels.reshape(image.shape[:2])
|
||||
@@ -0,0 +1,114 @@
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
import scipy
|
||||
cimport cython
|
||||
|
||||
from skimage.morphology.ccomp cimport find_root, join_trees
|
||||
|
||||
from ..util import img_as_float
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def _felzenszwalb_grey(image, double scale=1, sigma=0.8, int min_size=20):
|
||||
"""Felzenszwalb's efficient graph based segmentation for a single channel.
|
||||
|
||||
Produces an oversegmentation of a 2d image using a fast, minimum spanning
|
||||
tree based clustering on the image grid.
|
||||
The number of produced segments as well as their size can only be
|
||||
controlled indirectly through ``scale``. Segment size within an image can
|
||||
vary greatly depending on local contrast.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image: ndarray
|
||||
Input image.
|
||||
scale: float
|
||||
Sets the obervation level. Higher means larger clusters.
|
||||
sigma: float
|
||||
Width of Gaussian smoothing kernel used in preprocessing.
|
||||
Larger sigma gives smother segment boundaries.
|
||||
min_size: int
|
||||
Minimum component size. Enforced using postprocessing.
|
||||
|
||||
Returns
|
||||
-------
|
||||
segment_mask: (height, width) ndarray
|
||||
Integer mask indicating segment labels.
|
||||
"""
|
||||
if image.ndim != 2:
|
||||
raise ValueError("This algorithm works only on single-channel 2d"
|
||||
"images. Got image of shape %s" % str(image.shape))
|
||||
image = img_as_float(image)
|
||||
# rescale scale to behave like in reference implementation
|
||||
scale = float(scale) / 255.
|
||||
image = scipy.ndimage.gaussian_filter(image, sigma=sigma)
|
||||
|
||||
# compute edge weights in 8 connectivity:
|
||||
right_cost = np.abs((image[1:, :] - image[:-1, :]))
|
||||
down_cost = np.abs((image[:, 1:] - image[:, :-1]))
|
||||
dright_cost = np.abs((image[1:, 1:] - image[:-1, :-1]))
|
||||
uright_cost = np.abs((image[1:, :-1] - image[:-1, 1:]))
|
||||
cdef np.ndarray[np.float_t, ndim=1] costs = np.hstack([right_cost.ravel(),
|
||||
down_cost.ravel(), dright_cost.ravel(),
|
||||
uright_cost.ravel()]).astype(np.float)
|
||||
# compute edges between pixels:
|
||||
height, width = image.shape[:2]
|
||||
cdef np.ndarray[np.int_t, ndim=2] segments \
|
||||
= np.arange(width * height).reshape(height, width)
|
||||
right_edges = np.c_[segments[1:, :].ravel(), segments[:-1, :].ravel()]
|
||||
down_edges = np.c_[segments[:, 1:].ravel(), segments[:, :-1].ravel()]
|
||||
dright_edges = np.c_[segments[1:, 1:].ravel(), segments[:-1, :-1].ravel()]
|
||||
uright_edges = np.c_[segments[:-1, 1:].ravel(), segments[1:, :-1].ravel()]
|
||||
cdef np.ndarray[np.int_t, ndim=2] edges \
|
||||
= np.vstack([right_edges, down_edges, dright_edges, uright_edges])
|
||||
# initialize data structures for segment size
|
||||
# and inner cost, then start greedy iteration over edges.
|
||||
edge_queue = np.argsort(costs)
|
||||
edges = np.ascontiguousarray(edges[edge_queue])
|
||||
costs = np.ascontiguousarray(costs[edge_queue])
|
||||
cdef np.int_t *segments_p = <np.int_t*>segments.data
|
||||
cdef np.int_t *edges_p = <np.int_t*>edges.data
|
||||
cdef np.float_t *costs_p = <np.float_t*>costs.data
|
||||
cdef np.ndarray[np.int_t, ndim=1] segment_size \
|
||||
= np.ones(width * height, dtype=np.int)
|
||||
# inner cost of segments
|
||||
cdef np.ndarray[np.float_t, ndim=1] cint = np.zeros(width * height)
|
||||
cdef int seg0, seg1, seg_new, e
|
||||
cdef float cost, inner_cost0, inner_cost1
|
||||
# set costs_p back one. we increase it before we use it
|
||||
# since we might continue before that.
|
||||
costs_p -= 1
|
||||
for e in range(costs.size):
|
||||
seg0 = find_root(segments_p, edges_p[0])
|
||||
seg1 = find_root(segments_p, edges_p[1])
|
||||
edges_p += 2
|
||||
costs_p += 1
|
||||
if seg0 == seg1:
|
||||
continue
|
||||
inner_cost0 = cint[seg0] + scale / segment_size[seg0]
|
||||
inner_cost1 = cint[seg1] + scale / segment_size[seg1]
|
||||
if costs_p[0] < min(inner_cost0, inner_cost1):
|
||||
# update size and cost
|
||||
join_trees(segments_p, seg0, seg1)
|
||||
seg_new = find_root(segments_p, seg0)
|
||||
segment_size[seg_new] = segment_size[seg0] + segment_size[seg1]
|
||||
cint[seg_new] = costs_p[0]
|
||||
|
||||
# postprocessing to remove small segments
|
||||
edges_p = <np.int_t*>edges.data
|
||||
for e in range(costs.size):
|
||||
seg0 = find_root(segments_p, edges_p[0])
|
||||
seg1 = find_root(segments_p, edges_p[1])
|
||||
edges_p += 2
|
||||
if segment_size[seg0] < min_size or segment_size[seg1] < min_size:
|
||||
join_trees(segments_p, seg0, seg1)
|
||||
|
||||
# unravel the union find tree
|
||||
flat = segments.ravel()
|
||||
old = np.zeros_like(flat)
|
||||
while (old != flat).any():
|
||||
old = flat
|
||||
flat = flat[flat]
|
||||
flat = np.unique(flat, return_inverse=True)[1]
|
||||
return flat.reshape((height, width))
|
||||
@@ -0,0 +1,163 @@
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
cimport cython
|
||||
|
||||
from itertools import product
|
||||
from scipy import ndimage
|
||||
|
||||
from ..util import img_as_float
|
||||
from ..color import rgb2lab
|
||||
|
||||
|
||||
cdef extern from "math.h":
|
||||
double exp(double)
|
||||
double sqrt(double)
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def quickshift(image, ratio=1., float 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, between 0 and 1.
|
||||
Balances color-space proximity and image-space proximity.
|
||||
Higher values give more weight to color-space.
|
||||
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.
|
||||
sigma : float
|
||||
Width for Gaussian smoothing as preprocessing. Zero means no smoothing.
|
||||
convert2lab : bool
|
||||
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 or 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 = ndimage.gaussian_filter(img_as_float(image), [sigma, sigma, 0])
|
||||
cdef np.ndarray[dtype=np.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 int w = int(3 * kernel_size)
|
||||
|
||||
cdef int height = image_c.shape[0]
|
||||
cdef int width = image_c.shape[1]
|
||||
cdef int channels = image_c.shape[2]
|
||||
cdef double current_density, closest, dist
|
||||
|
||||
cdef int r, c, r_, c_, channel
|
||||
|
||||
cdef np.float_t* image_p = <np.float_t*> image_c.data
|
||||
cdef np.float_t* current_pixel_p = image_p
|
||||
|
||||
cdef np.ndarray[dtype=np.float_t, ndim=2] densities \
|
||||
= np.zeros((height, width))
|
||||
# compute densities
|
||||
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**2))
|
||||
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 np.ndarray[dtype=np.int_t, ndim=2] parent \
|
||||
= np.arange(width * height).reshape(height, width)
|
||||
cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent \
|
||||
= np.zeros((height, width))
|
||||
# find nearest node with higher density
|
||||
current_pixel_p = image_p
|
||||
for r in range(height):
|
||||
for c in range(width):
|
||||
current_density = densities[r, c]
|
||||
closest = np.inf
|
||||
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
|
||||
@@ -0,0 +1,124 @@
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
from time import time
|
||||
from scipy import ndimage
|
||||
from ..util import img_as_float
|
||||
from ..color import rgb2lab
|
||||
|
||||
|
||||
def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
|
||||
convert2lab=True):
|
||||
"""Segments image using k-means clustering in Color-(x,y) space.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : (width, height, 3) ndarray
|
||||
Input image.
|
||||
ratio: float
|
||||
Balances color-space proximity and image-space proximity.
|
||||
Higher values give more weight to color-space.
|
||||
max_iter : int
|
||||
Maximum number of iterations of k-means.
|
||||
sigma : float
|
||||
Width of Gaussian smoothing kernel for preprocessing. Zero means no
|
||||
smoothing.
|
||||
convert2lab : bool
|
||||
Whether the input should be converted to Lab colorspace prior to
|
||||
segmentation. For this purpose, the input is assumed to be RGB. Highly
|
||||
recommended.
|
||||
|
||||
Returns
|
||||
-------
|
||||
segment_mask : (width, height) ndarray
|
||||
Integer mask indicating segment labels.
|
||||
|
||||
Notes
|
||||
-----
|
||||
The image is smoothed using a Gaussian kernel prior to segmentation.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi,
|
||||
Pascal Fua, and Sabine Süsstrunk, SLIC Superpixels Compared to
|
||||
State-of-the-art Superpixel Methods, TPAMI, May 2012.
|
||||
|
||||
"""
|
||||
image = np.atleast_3d(image)
|
||||
if image.shape[2] != 3:
|
||||
ValueError("Only 3-channel 2D images are supported.")
|
||||
image = ndimage.gaussian_filter(img_as_float(image), [sigma, sigma, 0])
|
||||
if convert2lab:
|
||||
image = rgb2lab(image)
|
||||
|
||||
# initialize on grid:
|
||||
cdef int height, width
|
||||
height, width = image.shape[:2]
|
||||
# approximate grid size for desired n_segments
|
||||
cdef int step = np.ceil(np.sqrt(height * width / n_segments))
|
||||
grid_y, grid_x = np.mgrid[:height, :width]
|
||||
means_y = grid_y[::step, ::step]
|
||||
means_x = grid_x[::step, ::step]
|
||||
|
||||
means_color = np.zeros((means_y.shape[0], means_y.shape[1], 3))
|
||||
cdef np.ndarray[dtype=np.float_t, ndim=2] means \
|
||||
= np.dstack([means_y, means_x, means_color]).reshape(-1, 5)
|
||||
cdef np.float_t* current_mean
|
||||
cdef np.float_t* mean_entry
|
||||
n_means = means.shape[0]
|
||||
# we do the scaling of ratio in the same way as in the SLIC paper
|
||||
# so the values have the same meaning
|
||||
ratio = (ratio / float(step)) ** 2
|
||||
cdef np.ndarray[dtype=np.float_t, ndim=3] image_yx \
|
||||
= np.dstack([grid_y, grid_x, image / ratio]).copy("C")
|
||||
cdef int i, k, x, y, x_min, x_max, y_min, y_max, changes
|
||||
cdef double dist_mean
|
||||
|
||||
cdef np.ndarray[dtype=np.int_t, ndim=2] nearest_mean \
|
||||
= np.zeros((height, width), dtype=np.int)
|
||||
cdef np.ndarray[dtype=np.float_t, ndim=2] distance \
|
||||
= np.empty((height, width))
|
||||
cdef np.float_t* image_p = <np.float_t*> image_yx.data
|
||||
cdef np.float_t* distance_p = <np.float_t*> distance.data
|
||||
cdef np.float_t* current_distance
|
||||
cdef np.float_t* current_pixel
|
||||
cdef double tmp
|
||||
for i in range(max_iter):
|
||||
distance.fill(np.inf)
|
||||
changes = 0
|
||||
current_mean = <np.float_t*> means.data
|
||||
# assign pixels to means
|
||||
for k in range(n_means):
|
||||
# compute windows:
|
||||
y_min = int(max(current_mean[0] - 2 * step, 0))
|
||||
y_max = int(min(current_mean[0] + 2 * step, height))
|
||||
x_min = int(max(current_mean[1] - 2 * step, 0))
|
||||
x_max = int(min(current_mean[1] + 2 * step, width))
|
||||
for y in range(y_min, y_max):
|
||||
current_pixel = &image_p[5 * (y * width + x_min)]
|
||||
current_distance = &distance_p[y * width + x_min]
|
||||
for x in range(x_min, x_max):
|
||||
mean_entry = current_mean
|
||||
dist_mean = 0
|
||||
for c in range(5):
|
||||
# you would think the compiler can optimize the squaring
|
||||
# itself. mine can't (with O2)
|
||||
tmp = current_pixel[0] - mean_entry[0]
|
||||
dist_mean += tmp * tmp
|
||||
current_pixel += 1
|
||||
mean_entry += 1
|
||||
# some precision issue here. Doesnt work if testing ">"
|
||||
if current_distance[0] - dist_mean > 1e-10:
|
||||
nearest_mean[y, x] = k
|
||||
current_distance[0] = dist_mean
|
||||
changes += 1
|
||||
current_distance += 1
|
||||
current_mean += 5
|
||||
if changes == 0:
|
||||
break
|
||||
# recompute means:
|
||||
means_list = [np.bincount(nearest_mean.ravel(),
|
||||
image_yx[:, :, j].ravel()) for j in range(5)]
|
||||
in_mean = np.bincount(nearest_mean.ravel())
|
||||
in_mean[in_mean == 0] = 1
|
||||
means = (np.vstack(means_list) / in_mean).T.copy("C")
|
||||
return nearest_mean
|
||||
@@ -0,0 +1,19 @@
|
||||
import numpy as np
|
||||
from ..morphology import dilation, square
|
||||
from ..util import img_as_float
|
||||
|
||||
|
||||
def find_boundaries(label_img):
|
||||
boundaries = np.zeros(label_img.shape, dtype=np.bool)
|
||||
boundaries[1:, :] += label_img[1:, :] != label_img[:-1, :]
|
||||
boundaries[:, 1:] += label_img[:, 1:] != label_img[:, :-1]
|
||||
return boundaries
|
||||
|
||||
|
||||
def visualize_boundaries(img, label_img):
|
||||
img = img_as_float(img, force_copy=True)
|
||||
boundaries = find_boundaries(label_img)
|
||||
outer_boundaries = dilation(boundaries.astype(np.uint8), square(2))
|
||||
img[outer_boundaries != 0, :] = np.array([0, 0, 0]) # black
|
||||
img[boundaries, :] = np.array([1, 1, 0]) # yellow
|
||||
return img
|
||||
@@ -0,0 +1,34 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
from skimage._build import cython
|
||||
|
||||
base_path = os.path.abspath(os.path.dirname(__file__))
|
||||
|
||||
|
||||
def configuration(parent_package='', top_path=None):
|
||||
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
|
||||
|
||||
config = Configuration('segmentation', parent_package, top_path)
|
||||
|
||||
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'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
cython(['_slic.pyx'], working_path=base_path)
|
||||
config.add_extension('_slic', sources=['_slic.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.distutils.core import setup
|
||||
setup(maintainer='scikits-image Developers',
|
||||
maintainer_email='scikits-image@googlegroups.com',
|
||||
description='Segmentation Algorithms',
|
||||
url='https://github.com/scikits-image/scikits-image',
|
||||
license='SciPy License (BSD Style)',
|
||||
**(configuration(top_path='').todict())
|
||||
)
|
||||
@@ -0,0 +1,39 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_equal, assert_array_equal
|
||||
from nose.tools import assert_greater
|
||||
from skimage.segmentation import felzenszwalb
|
||||
|
||||
|
||||
def test_grey():
|
||||
# very weak tests. This algorithm is pretty unstable.
|
||||
img = np.zeros((20, 21))
|
||||
img[:10, 10:] = 0.2
|
||||
img[10:, :10] = 0.4
|
||||
img[10:, 10:] = 0.6
|
||||
seg = felzenszwalb(img, sigma=0)
|
||||
# we expect 4 segments:
|
||||
assert_equal(len(np.unique(seg)), 4)
|
||||
# that mostly respect the 4 regions:
|
||||
for i in xrange(4):
|
||||
hist = np.histogram(img[seg == i], bins=[0, 0.1, 0.3, 0.5, 1])[0]
|
||||
assert_greater(hist[i], 40)
|
||||
|
||||
|
||||
def test_color():
|
||||
# very weak tests. This algorithm is pretty unstable.
|
||||
img = np.zeros((20, 21, 3))
|
||||
img[:10, :10, 0] = 1
|
||||
img[10:, :10, 1] = 1
|
||||
img[10:, 10:, 2] = 1
|
||||
seg = felzenszwalb(img, sigma=0)
|
||||
# we expect 4 segments:
|
||||
assert_equal(len(np.unique(seg)), 4)
|
||||
assert_array_equal(seg[:10, :10], 0)
|
||||
assert_array_equal(seg[10:, :10], 2)
|
||||
assert_array_equal(seg[:10, 10:], 1)
|
||||
assert_array_equal(seg[10:, 10:], 3)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
@@ -0,0 +1,52 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_equal, assert_array_equal
|
||||
from nose.tools import assert_true, assert_greater
|
||||
from skimage.segmentation import quickshift
|
||||
|
||||
|
||||
def test_grey():
|
||||
rnd = np.random.RandomState(0)
|
||||
img = np.zeros((20, 21))
|
||||
img[:10, 10:] = 0.2
|
||||
img[10:, :10] = 0.4
|
||||
img[10:, 10:] = 0.6
|
||||
img += 0.1 * rnd.normal(size=img.shape)
|
||||
seg = quickshift(img, kernel_size=2, max_dist=3, random_seed=0,
|
||||
convert2lab=False, sigma=0)
|
||||
# we expect 4 segments:
|
||||
assert_equal(len(np.unique(seg)), 4)
|
||||
# that mostly respect the 4 regions:
|
||||
for i in xrange(4):
|
||||
hist = np.histogram(img[seg == i], bins=[0, 0.1, 0.3, 0.5, 1])[0]
|
||||
assert_greater(hist[i], 20)
|
||||
|
||||
|
||||
def test_color():
|
||||
rnd = np.random.RandomState(0)
|
||||
img = np.zeros((20, 21, 3))
|
||||
img[:10, :10, 0] = 1
|
||||
img[10:, :10, 1] = 1
|
||||
img[10:, 10:, 2] = 1
|
||||
img += 0.01 * rnd.normal(size=img.shape)
|
||||
img[img > 1] = 1
|
||||
img[img < 0] = 0
|
||||
seg = quickshift(img, random_seed=0, max_dist=30, kernel_size=10, sigma=0)
|
||||
# we expect 4 segments:
|
||||
assert_equal(len(np.unique(seg)), 4)
|
||||
assert_array_equal(seg[:10, :10], 0)
|
||||
assert_array_equal(seg[10:, :10], 3)
|
||||
assert_array_equal(seg[:10, 10:], 1)
|
||||
assert_array_equal(seg[10:, 10:], 2)
|
||||
|
||||
seg2 = quickshift(img, kernel_size=1, max_dist=2, random_seed=0,
|
||||
convert2lab=False, sigma=0)
|
||||
# very oversegmented:
|
||||
assert_equal(len(np.unique(seg2)), 7)
|
||||
# still don't cross lines
|
||||
assert_true((seg2[9, :] != seg2[10, :]).all())
|
||||
assert_true((seg2[:, 9] != seg2[:, 10]).all())
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
@@ -0,0 +1,27 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_equal, assert_array_equal
|
||||
from skimage.segmentation import slic
|
||||
|
||||
|
||||
def test_color():
|
||||
rnd = np.random.RandomState(0)
|
||||
img = np.zeros((20, 21, 3))
|
||||
img[:10, :10, 0] = 1
|
||||
img[10:, :10, 1] = 1
|
||||
img[10:, 10:, 2] = 1
|
||||
img += 0.01 * rnd.normal(size=img.shape)
|
||||
img[img > 1] = 1
|
||||
img[img < 0] = 0
|
||||
seg = slic(img, sigma=0, n_segments=4)
|
||||
# we expect 4 segments:
|
||||
print(seg)
|
||||
assert_equal(len(np.unique(seg)), 4)
|
||||
assert_array_equal(seg[:10, :10], 0)
|
||||
assert_array_equal(seg[10:, :10], 2)
|
||||
assert_array_equal(seg[:10, 10:], 1)
|
||||
assert_array_equal(seg[10:, 10:], 3)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
@@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_subpackage('morphology')
|
||||
config.add_subpackage('transform')
|
||||
config.add_subpackage('util')
|
||||
config.add_subpackage('segmentation')
|
||||
|
||||
def add_test_directories(arg, dirname, fnames):
|
||||
if dirname.split(os.path.sep)[-1] == 'tests':
|
||||
|
||||
Reference in New Issue
Block a user