From b8c0663332398d0487c3e43d205cb02b478b243d Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 15 Jun 2012 20:59:17 +0200 Subject: [PATCH 01/63] Add segmentation setup.py for felsenzwalb algorithm --- skimage/segmentation/felsenzwalb.pyx | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 skimage/segmentation/felsenzwalb.pyx diff --git a/skimage/segmentation/felsenzwalb.pyx b/skimage/segmentation/felsenzwalb.pyx new file mode 100644 index 00000000..38a95411 --- /dev/null +++ b/skimage/segmentation/felsenzwalb.pyx @@ -0,0 +1,4 @@ +# Implements Felsenzwalb's efficient graph based image segmentation. +# Author: Andreas Mueller + +def felsenzwalb(np.ndarray[dtype= From 967eb5b50d959543bc41904fbdf651418acbac17 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 16:42:18 +0200 Subject: [PATCH 02/63] ENH first draft of felzenszwalbs graph based image segmentation in Python --- .../plot_felzenszwalb_segmentation.py | 13 ++++ skimage/segmentation/__init__.py | 1 + skimage/segmentation/felzenszwalb.py | 60 +++++++++++++++++++ 3 files changed, 74 insertions(+) create mode 100644 doc/examples/plot_felzenszwalb_segmentation.py create mode 100644 skimage/segmentation/felzenszwalb.py diff --git a/doc/examples/plot_felzenszwalb_segmentation.py b/doc/examples/plot_felzenszwalb_segmentation.py new file mode 100644 index 00000000..dbcb2abb --- /dev/null +++ b/doc/examples/plot_felzenszwalb_segmentation.py @@ -0,0 +1,13 @@ +import matplotlib.pyplot as plt +import numpy as np + +from skimage.data import lena +from skimage.segmentation import felzenszwalb_segmentation + +img = lena() +segments = felzenszwalb_segmentation(img, k=1000) +plt.imshow(img) +plt.figure() +plt.imshow(segments) +plt.show() +print("num segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 118c3ec2..372c58fc 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1 +1,2 @@ from .random_walker_segmentation import random_walker +from .felzenszwalb import felzenszwalb_segmentation diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py new file mode 100644 index 00000000..9dea47a6 --- /dev/null +++ b/skimage/segmentation/felzenszwalb.py @@ -0,0 +1,60 @@ +import numpy as np +from collections import defaultdict +import scipy + +#from ..util import img_as_float +#from ..color import rgb2grey +from .union_find import UnionFind + +from IPython.core.debugger import Tracer +tracer = Tracer() + + +def felzenszwalb_segmentation(image, k, sigma=0.8): + k = float(k) + #image = img_as_float(image) + #image = rgb2grey(image) + image = image[:, :, 0] + image = scipy.ndimage.gaussian_filter(image, sigma=sigma) + + # compute edge weights in 8 connectivity: + #right_cost = np.sum((image[1:, :, :] - image[:-1, :, :]) ** 2, axis=2) + #down_cost = np.sum((image[:, 1:, :] - image[:, :-1, :]) ** 2, axis=2) + 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:])) + costs = np.hstack([right_cost.ravel(), down_cost.ravel(), + dright_cost.ravel(), uright_cost.ravel()]) + # compute edges between pixels: + width, height = image.shape[:2] + indices = np.arange(width * height).reshape(width, height) + right_edges = np.c_[indices[1:, :].ravel(), indices[:-1, :].ravel()] + down_edges = np.c_[indices[:, 1:].ravel(), indices[:, :-1].ravel()] + dright_edges = np.c_[indices[1:, 1:].ravel(), indices[:-1, :-1].ravel()] + uright_edges = np.c_[indices[:-1, 1:].ravel(), indices[1:, :-1].ravel()] + 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) + segments = UnionFind() + segment_size = defaultdict(lambda: 1) + # inner cost of segments + cint = defaultdict(lambda: 0) + for edge, cost in zip(edges[edge_queue], costs[edge_queue]): + seg0 = segments[edge[0]] + seg1 = segments[edge[1]] + if seg0 == seg1: + continue + inner_cost0 = cint[seg0] + k / segment_size[seg0] + inner_cost1 = cint[seg1] + k / segment_size[seg1] + if cost < min(inner_cost0, inner_cost1): + seg_new = segments.union(seg0, seg1) + # update size and cost + segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] + cint[seg_new] = cost + out = np.zeros(width * height, dtype=np.int) + for i in xrange(width * height): + out[i] = segments[i] + out = out.reshape(width, height) + return out From e2d60f01357d9d18ee3c34be4b7953e7dd634168 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 17:28:43 +0200 Subject: [PATCH 03/63] ENH using union find from morphology module --- skimage/morphology/ccomp.pxd | 10 +++++++ skimage/segmentation/felsenzwalb.pyx | 4 --- .../{felzenszwalb.py => felzenszwalb.pyx} | 20 +++++++------- skimage/segmentation/setup.py | 27 +++++++++++++++++++ skimage/setup.py | 1 + 5 files changed, 48 insertions(+), 14 deletions(-) create mode 100644 skimage/morphology/ccomp.pxd delete mode 100644 skimage/segmentation/felsenzwalb.pyx rename skimage/segmentation/{felzenszwalb.py => felzenszwalb.pyx} (82%) create mode 100644 skimage/segmentation/setup.py diff --git a/skimage/morphology/ccomp.pxd b/skimage/morphology/ccomp.pxd new file mode 100644 index 00000000..0b431832 --- /dev/null +++ b/skimage/morphology/ccomp.pxd @@ -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) diff --git a/skimage/segmentation/felsenzwalb.pyx b/skimage/segmentation/felsenzwalb.pyx deleted file mode 100644 index 38a95411..00000000 --- a/skimage/segmentation/felsenzwalb.pyx +++ /dev/null @@ -1,4 +0,0 @@ -# Implements Felsenzwalb's efficient graph based image segmentation. -# Author: Andreas Mueller - -def felsenzwalb(np.ndarray[dtype= diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.pyx similarity index 82% rename from skimage/segmentation/felzenszwalb.py rename to skimage/segmentation/felzenszwalb.pyx index 9dea47a6..5a0ecbac 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.pyx @@ -1,10 +1,11 @@ import numpy as np +cimport numpy as np from collections import defaultdict import scipy #from ..util import img_as_float #from ..color import rgb2grey -from .union_find import UnionFind +from skimage.morphology.ccomp cimport find_root, join_trees from IPython.core.debugger import Tracer tracer = Tracer() @@ -37,24 +38,23 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): # initialize data structures for segment size # and inner cost, then start greedy iteration over edges. edge_queue = np.argsort(costs) - segments = UnionFind() + cdef np.ndarray[np.int_t, ndim=2] segments = indices.reshape(width, height) + cdef np.int_t *segments_p = segments.data + cdef np.int_t seg_new segment_size = defaultdict(lambda: 1) # inner cost of segments cint = defaultdict(lambda: 0) for edge, cost in zip(edges[edge_queue], costs[edge_queue]): - seg0 = segments[edge[0]] - seg1 = segments[edge[1]] + seg0 = find_root(segments_p, edge[0]) + seg1 = find_root(segments_p, edge[1]) if seg0 == seg1: continue inner_cost0 = cint[seg0] + k / segment_size[seg0] inner_cost1 = cint[seg1] + k / segment_size[seg1] if cost < min(inner_cost0, inner_cost1): - seg_new = segments.union(seg0, seg1) # 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] = cost - out = np.zeros(width * height, dtype=np.int) - for i in xrange(width * height): - out[i] = segments[i] - out = out.reshape(width, height) - return out + return segments diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py new file mode 100644 index 00000000..8f6361bf --- /dev/null +++ b/skimage/segmentation/setup.py @@ -0,0 +1,27 @@ +#!/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.pyx'], working_path=base_path) + config.add_extension('felzenszwalb', sources=['felzenszwalb.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()) + ) diff --git a/skimage/setup.py b/skimage/setup.py index 0afc5bc9..02c0b52d 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -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': From b1b1c343b452b32371db7a5c6422afd8fce16c4c Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 19:15:20 +0200 Subject: [PATCH 04/63] MISC remove debugging tracer, unnecessary variable. --- skimage/segmentation/felzenszwalb.pyx | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/felzenszwalb.pyx index 5a0ecbac..a53bf70e 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/felzenszwalb.pyx @@ -7,9 +7,6 @@ import scipy #from ..color import rgb2grey from skimage.morphology.ccomp cimport find_root, join_trees -from IPython.core.debugger import Tracer -tracer = Tracer() - def felzenszwalb_segmentation(image, k, sigma=0.8): k = float(k) @@ -29,16 +26,15 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): dright_cost.ravel(), uright_cost.ravel()]) # compute edges between pixels: width, height = image.shape[:2] - indices = np.arange(width * height).reshape(width, height) - right_edges = np.c_[indices[1:, :].ravel(), indices[:-1, :].ravel()] - down_edges = np.c_[indices[:, 1:].ravel(), indices[:, :-1].ravel()] - dright_edges = np.c_[indices[1:, 1:].ravel(), indices[:-1, :-1].ravel()] - uright_edges = np.c_[indices[:-1, 1:].ravel(), indices[1:, :-1].ravel()] + cdef np.ndarray[np.int_t, ndim=2] segments = np.arange(width * height).reshape(width, height) + 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()] 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) - cdef np.ndarray[np.int_t, ndim=2] segments = indices.reshape(width, height) cdef np.int_t *segments_p = segments.data cdef np.int_t seg_new segment_size = defaultdict(lambda: 1) From 40ecdd29dbce103f563de7a62f6628a45a3d61ac Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 19:15:43 +0200 Subject: [PATCH 05/63] ENH naive pure python implementation of quickshift --- skimage/segmentation/quickshift.py | 47 ++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 skimage/segmentation/quickshift.py diff --git a/skimage/segmentation/quickshift.py b/skimage/segmentation/quickshift.py new file mode 100644 index 00000000..5025b3a2 --- /dev/null +++ b/skimage/segmentation/quickshift.py @@ -0,0 +1,47 @@ +import numpy as np +from itertools import product, combinations_with_replacement + +from IPython.core.debugger import Tracer +tracer = Tracer() + + +def quickshift(image, sigma=5, tau=10): + # do smoothing beforehand? + width, height = image.shape[:2] + densities = np.zeros((width, height)) + w = 10 + + # TODO: normalize density by number of considered points. + # important for the border! + # compute densities + for x, y in product(xrange(width), xrange(height)): + current_pixel = np.hstack([image[x, y, :], x, y]) + for xx, yy in combinations_with_replacement(xrange(-w / 2, w / 2), 2): + x_, y_ = x + xx, y + yy + if 0 <= x_ < width and 0 <= y_ < height: + other_pixel = np.hstack([image[x_, y_, :], x_, y_]) + dist = np.sum((current_pixel - other_pixel) ** 2) + densities[x, y] += np.exp(-dist / sigma) + + # default parent to self: + parent = np.arange(width * height).reshape(width, height) + # find nearest node with higher density + for x, y in product(xrange(width), xrange(height)): + current_density = densities[x, y] + current_pixel = np.hstack([image[x, y, :], x, y]) + closest = np.inf + for xx, yy in combinations_with_replacement(xrange(-w / 2, w / 2), 2): + x_, y_ = x + xx, y + yy + if 0 <= x_ < width and 0 <= y_ < height: + if densities[x_, y_] > current_density: + other_pixel = np.hstack([image[x_, y_, :], x_, y_]) + dist = np.sum((current_pixel - other_pixel) ** 2) + if dist < closest: + closest = dist + parent[x, y] = x_ * width + y_ + flat = parent.ravel() + old = np.zeros_like(flat) + while (old != flat).any(): + old = flat + flat = flat[flat] + return flat.reshape(parent.shape) From eb5c2fe5d49ccf93c9499c000c17e6a8fbb45bbf Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 20:52:04 +0200 Subject: [PATCH 06/63] ENH fixed stupid bug in quickshift, example --- doc/examples/plot_quickshift.py | 47 ++++++++++++++++++++++++++++ skimage/segmentation/__init__.py | 5 ++- skimage/segmentation/quickshift.py | 50 ++++++++++++++++++++++++------ 3 files changed, 91 insertions(+), 11 deletions(-) create mode 100644 doc/examples/plot_quickshift.py diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py new file mode 100644 index 00000000..80f216fa --- /dev/null +++ b/doc/examples/plot_quickshift.py @@ -0,0 +1,47 @@ +import matplotlib.pyplot as plt +import numpy as np + +from scipy import ndimage +#from skimage.data import lena +#from skimage.util import img_as_float +from skimage.segmentation import quickshift + +from IPython.core.debugger import Tracer +tracer = Tracer() + + +def microstructure(l=256): + """ + Synthetic binary data: binary microstructure with blobs. + + Parameters + ---------- + + l: int, optional + linear size of the returned image + """ + n = 5 + x, y = np.ogrid[0:l, 0:l] + mask = np.zeros((l, l)) + generator = np.random.RandomState(1) + points = l * generator.rand(2, n ** 2) + mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 + mask = ndimage.gaussian_filter(mask, sigma=l / (4. * n)) + return (mask > mask.mean()).astype(np.float) + + +#img = img_as_float(lena()[250:300, 250:300]) +img = microstructure(l=50) +segments = quickshift(img.reshape(50, 50, 1)) +segments = np.unique(segments, return_inverse=True)[1].reshape(50, 50) +intensities = np.bincount(segments.ravel(), img.ravel()) +counts = np.bincount(segments.ravel()) +intensities /= counts + +plt.imshow(img, interpolation='nearest') +plt.figure() +plt.imshow(segments, interpolation='nearest') +plt.figure() +plt.imshow(intensities[segments], interpolation='nearest') +plt.show() +print("num segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 372c58fc..0ea91444 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,2 +1,5 @@ from .random_walker_segmentation import random_walker -from .felzenszwalb import felzenszwalb_segmentation +#from .felzenszwalb import felzenszwalb_segmentation +from .quickshift import quickshift + +__all__ = [random_walker, quickshift] diff --git a/skimage/segmentation/quickshift.py b/skimage/segmentation/quickshift.py index 5025b3a2..65c8847f 100644 --- a/skimage/segmentation/quickshift.py +++ b/skimage/segmentation/quickshift.py @@ -1,36 +1,62 @@ import numpy as np -from itertools import product, combinations_with_replacement - -from IPython.core.debugger import Tracer -tracer = Tracer() +from itertools import product def quickshift(image, sigma=5, tau=10): - # do smoothing beforehand? + """Computes quickshift clustering in RGB-(x,y) space. + + Parameters + ---------- + image: ndarray, [width, height, channels] + Input image + sigma: float + Width of Gaussian kernel used in smoothing the + sample density. Higher means less clusters. + tau: float + Cut-off point for data distances. + Higher means less clusters. + + Returns + ------- + segment_mask: ndarray, [width, height] + Integer mask indicating segment labels. + """ + + # We compute the distances twice since otherwise + # we might get crazy memory overhead (width * height * windowsize**2) + + # TODO do smoothing beforehand? + # TODO manage borders somehow? + + # window size for neighboring pixels to consider + if sigma < 1: + raise ValueError("Sigma should be >= 1") + w = int(2 * sigma) + width, height = image.shape[:2] densities = np.zeros((width, height)) - w = 10 - # TODO: normalize density by number of considered points. - # important for the border! # compute densities for x, y in product(xrange(width), xrange(height)): current_pixel = np.hstack([image[x, y, :], x, y]) - for xx, yy in combinations_with_replacement(xrange(-w / 2, w / 2), 2): + for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: other_pixel = np.hstack([image[x_, y_, :], x_, y_]) dist = np.sum((current_pixel - other_pixel) ** 2) densities[x, y] += np.exp(-dist / sigma) + # this will break ties that otherwise would give us headache + densities += np.random.normal(scale=0.00001, size=densities.shape) # default parent to self: parent = np.arange(width * height).reshape(width, height) + dist_parent = np.zeros((width, height)) # find nearest node with higher density for x, y in product(xrange(width), xrange(height)): current_density = densities[x, y] current_pixel = np.hstack([image[x, y, :], x, y]) closest = np.inf - for xx, yy in combinations_with_replacement(xrange(-w / 2, w / 2), 2): + for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: if densities[x_, y_] > current_density: @@ -39,7 +65,11 @@ def quickshift(image, sigma=5, tau=10): if dist < closest: closest = dist parent[x, y] = x_ * width + y_ + dist_parent[x, y] = closest + + dist_parent = dist_parent.ravel() flat = parent.ravel() + flat[dist_parent > tau] = np.arange(width * height)[dist_parent > tau] old = np.zeros_like(flat) while (old != flat).any(): old = flat From 8c735b64708c0d73356d096e51a536c888124596 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 21:09:07 +0200 Subject: [PATCH 07/63] ENH start cythonizing quickshift, get rid of hstack. --- .../{quickshift.py => quickshift.pyx} | 25 +++++++++++-------- skimage/segmentation/setup.py | 18 +++++++------ 2 files changed, 25 insertions(+), 18 deletions(-) rename skimage/segmentation/{quickshift.py => quickshift.pyx} (77%) diff --git a/skimage/segmentation/quickshift.py b/skimage/segmentation/quickshift.pyx similarity index 77% rename from skimage/segmentation/quickshift.py rename to skimage/segmentation/quickshift.pyx index 65c8847f..8bec6d2b 100644 --- a/skimage/segmentation/quickshift.py +++ b/skimage/segmentation/quickshift.pyx @@ -1,8 +1,10 @@ import numpy as np +cimport numpy as np + from itertools import product -def quickshift(image, sigma=5, tau=10): +def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, tau=10): """Computes quickshift clustering in RGB-(x,y) space. Parameters @@ -31,37 +33,38 @@ def quickshift(image, sigma=5, tau=10): # window size for neighboring pixels to consider if sigma < 1: raise ValueError("Sigma should be >= 1") - w = int(2 * sigma) + cdef int w = int(2 * sigma) + + cdef int width = image.shape[0] + cdef int height = image.shape[1] - width, height = image.shape[:2] - densities = np.zeros((width, height)) + cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) # compute densities for x, y in product(xrange(width), xrange(height)): - current_pixel = np.hstack([image[x, y, :], x, y]) + current_pixel = image[x, y, :] for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: - other_pixel = np.hstack([image[x_, y_, :], x_, y_]) - dist = np.sum((current_pixel - other_pixel) ** 2) + dist = np.sum((current_pixel - image[x_, y_, :])**2) + (x - x_)**2 + (y - y_)**2 densities[x, y] += np.exp(-dist / sigma) # this will break ties that otherwise would give us headache - densities += np.random.normal(scale=0.00001, size=densities.shape) + + densities += np.random.normal(scale=0.00001, size=(width, height)) # default parent to self: parent = np.arange(width * height).reshape(width, height) dist_parent = np.zeros((width, height)) # find nearest node with higher density for x, y in product(xrange(width), xrange(height)): current_density = densities[x, y] - current_pixel = np.hstack([image[x, y, :], x, y]) + current_pixel = image[x, y, :] closest = np.inf for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: if densities[x_, y_] > current_density: - other_pixel = np.hstack([image[x_, y_, :], x_, y_]) - dist = np.sum((current_pixel - other_pixel) ** 2) + dist = np.sum((current_pixel - image[x_, y_, :])**2) + (x - x_)**2 + (y - y_)**2 if dist < closest: closest = dist parent[x, y] = x_ * width + y_ diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index 8f6361bf..b6f8d1e2 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -5,23 +5,27 @@ 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.pyx'], working_path=base_path) - config.add_extension('felzenszwalb', sources=['felzenszwalb.c'], + #cython(['felzenszwalb.pyx'], working_path=base_path) + #config.add_extension('felzenszwalb', sources=['felzenszwalb.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()]) 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)', + 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()) ) From de52692a9d46fcd875869dff1c032b2ce7e0ce15 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 21:35:15 +0200 Subject: [PATCH 08/63] misc Uncomment Felzenszwalb as it is not messing with quickshift. --- skimage/segmentation/setup.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index b6f8d1e2..a197555f 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -11,9 +11,9 @@ def configuration(parent_package='', top_path=None): config = Configuration('segmentation', parent_package, top_path) - #cython(['felzenszwalb.pyx'], working_path=base_path) - #config.add_extension('felzenszwalb', sources=['felzenszwalb.c'], - #include_dirs=[get_numpy_include_dirs()]) + cython(['felzenszwalb.pyx'], working_path=base_path) + config.add_extension('felzenszwalb', sources=['felzenszwalb.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()]) From 48fa3252beb34ce31601830d8ad742cb334d4c71 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 21:53:28 +0200 Subject: [PATCH 09/63] ENH reasonable speed. --- skimage/segmentation/quickshift.pyx | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 8bec6d2b..5fd3eafe 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -3,6 +3,9 @@ cimport numpy as np from itertools import product +cdef extern from "math.h": + double exp(double) + def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, tau=10): """Computes quickshift clustering in RGB-(x,y) space. @@ -37,6 +40,9 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta cdef int width = image.shape[0] cdef int height = image.shape[1] + cdef int channels = image.shape[2] + cdef float closest, dist + cdef int x, y, xx, yy, x_, y_ cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) @@ -46,15 +52,18 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: - dist = np.sum((current_pixel - image[x_, y_, :])**2) + (x - x_)**2 + (y - y_)**2 - densities[x, y] += np.exp(-dist / sigma) + dist = 0 + for c in xrange(channels): + dist += (current_pixel[c] - image[x_, y_, c])**2 + dist += (x - x_)**2 + (y - y_)**2 + densities[x, y] += float(exp(-dist / sigma)) # this will break ties that otherwise would give us headache densities += np.random.normal(scale=0.00001, size=(width, height)) # default parent to self: - parent = np.arange(width * height).reshape(width, height) - dist_parent = np.zeros((width, height)) + cdef np.ndarray[dtype=np.int_t, ndim=2] parent = np.arange(width * height).reshape(width, height) + cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) # find nearest node with higher density for x, y in product(xrange(width), xrange(height)): current_density = densities[x, y] @@ -64,17 +73,20 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: if densities[x_, y_] > current_density: - dist = np.sum((current_pixel - image[x_, y_, :])**2) + (x - x_)**2 + (y - y_)**2 + dist = 0 + for c in xrange(channels): + dist += (current_pixel[c] - image[x_, y_, c])**2 + dist += (x - x_)**2 + (y - y_)**2 if dist < closest: closest = dist parent[x, y] = x_ * width + y_ dist_parent[x, y] = closest - dist_parent = dist_parent.ravel() + dist_parent_flat = dist_parent.ravel() flat = parent.ravel() - flat[dist_parent > tau] = np.arange(width * height)[dist_parent > tau] + flat[dist_parent_flat > tau] = np.arange(width * height)[dist_parent_flat > tau] old = np.zeros_like(flat) while (old != flat).any(): old = flat flat = flat[flat] - return flat.reshape(parent.shape) + return flat.reshape(width, height) From b977d59c1bd2cae4f91c49e23421f64c88972261 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 23:21:44 +0200 Subject: [PATCH 10/63] Color example :) --- doc/examples/plot_quickshift.py | 52 ++++++++++------------------- skimage/segmentation/quickshift.pyx | 18 ++++++++-- 2 files changed, 33 insertions(+), 37 deletions(-) diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index 80f216fa..08757e25 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -1,47 +1,31 @@ import matplotlib.pyplot as plt import numpy as np -from scipy import ndimage -#from skimage.data import lena -#from skimage.util import img_as_float +from skimage.data import lena from skimage.segmentation import quickshift +from skimage.util import img_as_float from IPython.core.debugger import Tracer tracer = Tracer() -def microstructure(l=256): - """ - Synthetic binary data: binary microstructure with blobs. - - Parameters - ---------- - - l: int, optional - linear size of the returned image - """ - n = 5 - x, y = np.ogrid[0:l, 0:l] - mask = np.zeros((l, l)) - generator = np.random.RandomState(1) - points = l * generator.rand(2, n ** 2) - mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 - mask = ndimage.gaussian_filter(mask, sigma=l / (4. * n)) - return (mask > mask.mean()).astype(np.float) - - -#img = img_as_float(lena()[250:300, 250:300]) -img = microstructure(l=50) -segments = quickshift(img.reshape(50, 50, 1)) -segments = np.unique(segments, return_inverse=True)[1].reshape(50, 50) -intensities = np.bincount(segments.ravel(), img.ravel()) -counts = np.bincount(segments.ravel()) -intensities /= counts +img = img_as_float(lena())[::3, ::3, :].copy("C") +segments = quickshift(img, sigma=2) +segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) +plt.subplot(131, title="original") plt.imshow(img, interpolation='nearest') -plt.figure() -plt.imshow(segments, interpolation='nearest') -plt.figure() -plt.imshow(intensities[segments], interpolation='nearest') + +plt.subplot(132, title="superpixels") +# shuffle the labels for better visualization +permuted_labels = np.random.permutation(segments.max() + 1) +plt.imshow(permuted_labels[segments], interpolation='nearest') + +plt.subplot(133, title="mean color") +colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in + xrange(img.shape[2])] +counts = np.bincount(segments.ravel()) +colors = np.vstack(colors) / counts +plt.imshow(colors.T[segments], interpolation='nearest') plt.show() print("num segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 5fd3eafe..3fbfbd98 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -3,11 +3,13 @@ cimport numpy as np from itertools import product +from time import time + cdef extern from "math.h": double exp(double) -def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, tau=10): +def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, tau=10, return_tree=False): """Computes quickshift clustering in RGB-(x,y) space. Parameters @@ -20,6 +22,8 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta tau: float Cut-off point for data distances. Higher means less clusters. + return_tree: bool + Whether to return the full segmentation hierarchy tree Returns ------- @@ -45,7 +49,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta cdef int x, y, xx, yy, x_, y_ cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) - + start = time() # compute densities for x, y in product(xrange(width), xrange(height)): current_pixel = image[x, y, :] @@ -57,6 +61,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta dist += (current_pixel[c] - image[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 densities[x, y] += float(exp(-dist / sigma)) + print("densities: %f" % (time() - start)) # this will break ties that otherwise would give us headache @@ -64,6 +69,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta # default parent to self: cdef np.ndarray[dtype=np.int_t, ndim=2] parent = np.arange(width * height).reshape(width, height) cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) + start = time() # find nearest node with higher density for x, y in product(xrange(width), xrange(height)): current_density = densities[x, y] @@ -81,7 +87,9 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta closest = dist parent[x, y] = x_ * width + y_ dist_parent[x, y] = closest + print("parents: %f" % (time() - start)) + start = time() dist_parent_flat = dist_parent.ravel() flat = parent.ravel() flat[dist_parent_flat > tau] = np.arange(width * height)[dist_parent_flat > tau] @@ -89,4 +97,8 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta while (old != flat).any(): old = flat flat = flat[flat] - return flat.reshape(width, height) + print("rest: %f" % (time() - start)) + flat = flat.reshape(width, height) + if return_tree: + return flat, parent + return flat From be4b44bc63a14b554ed80875ea00d2e8cae97666 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 16 Jun 2012 23:39:48 +0200 Subject: [PATCH 11/63] ENH CRAZY speedup --- doc/examples/plot_quickshift.py | 4 ++-- skimage/segmentation/quickshift.pyx | 15 +++++++++++---- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index 08757e25..59850469 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -9,8 +9,8 @@ from IPython.core.debugger import Tracer tracer = Tracer() -img = img_as_float(lena())[::3, ::3, :].copy("C") -segments = quickshift(img, sigma=2) +img = img_as_float(lena())[::2, ::2, :].copy("C") +segments = quickshift(img) segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) plt.subplot(131, title="original") diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 3fbfbd98..4a619fc2 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -47,20 +47,26 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta cdef int channels = image.shape[2] cdef float closest, dist cdef int x, y, xx, yy, x_, y_ + + cdef np.float_t* image_p = image.data + cdef np.float_t* current_pixel_p = image_p + cdef np.float_t* current_entry_p cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) start = time() # compute densities for x, y in product(xrange(width), xrange(height)): - current_pixel = image[x, y, :] for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy if 0 <= x_ < width and 0 <= y_ < height: dist = 0 + current_entry_p = current_pixel_p for c in xrange(channels): - dist += (current_pixel[c] - image[x_, y_, c])**2 + dist += (current_pixel_p[c] - image[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 densities[x, y] += float(exp(-dist / sigma)) + current_pixel_p += channels + print("densities: %f" % (time() - start)) # this will break ties that otherwise would give us headache @@ -71,9 +77,9 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) start = time() # find nearest node with higher density + current_pixel_p = image_p for x, y in product(xrange(width), xrange(height)): current_density = densities[x, y] - current_pixel = image[x, y, :] closest = np.inf for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): x_, y_ = x + xx, y + yy @@ -81,12 +87,13 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta if densities[x_, y_] > current_density: dist = 0 for c in xrange(channels): - dist += (current_pixel[c] - image[x_, y_, c])**2 + dist += (current_pixel_p[c] - image[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 if dist < closest: closest = dist parent[x, y] = x_ * width + y_ dist_parent[x, y] = closest + current_pixel_p += channels print("parents: %f" % (time() - start)) start = time() From cb3dba7847f8e01ed3de2e9b94082e8ad944c8c9 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 00:19:35 +0200 Subject: [PATCH 12/63] Bigger example --- doc/examples/plot_quickshift.py | 4 ++-- skimage/segmentation/quickshift.pyx | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index 59850469..b21c56e8 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -9,8 +9,8 @@ from IPython.core.debugger import Tracer tracer = Tracer() -img = img_as_float(lena())[::2, ::2, :].copy("C") -segments = quickshift(img) +img = img_as_float(lena()) +segments = quickshift(img, sigma=5, tau=20) segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) plt.subplot(131, title="original") diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 4a619fc2..7e27fc56 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -33,6 +33,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta # We compute the distances twice since otherwise # we might get crazy memory overhead (width * height * windowsize**2) + # if you want to speed up things: computing exp in C is the bottleneck ;) # TODO do smoothing beforehand? # TODO manage borders somehow? @@ -64,7 +65,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta for c in xrange(channels): dist += (current_pixel_p[c] - image[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 - densities[x, y] += float(exp(-dist / sigma)) + densities[x, y] += exp(-dist / sigma) current_pixel_p += channels print("densities: %f" % (time() - start)) From 58237a558a457ef7a32c67b99c2f81aa2b04e89a Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 20:52:45 +0200 Subject: [PATCH 13/63] ENH dirty fix, works though. Starting profiling. --- skimage/segmentation/__init__.py | 4 +-- skimage/segmentation/felzenszwalb.pyx | 52 +++++++++++++++++++++++++-- 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 0ea91444..8fa8dfb8 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,5 +1,5 @@ from .random_walker_segmentation import random_walker -#from .felzenszwalb import felzenszwalb_segmentation +from .felzenszwalb import felzenszwalb_segmentation from .quickshift import quickshift -__all__ = [random_walker, quickshift] +__all__ = [random_walker, quickshift, felzenszwalb_segmentation] diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/felzenszwalb.pyx index a53bf70e..11f3e93b 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/felzenszwalb.pyx @@ -5,7 +5,49 @@ import scipy #from ..util import img_as_float #from ..color import rgb2grey -from skimage.morphology.ccomp cimport find_root, join_trees +#from skimage.morphology.ccomp cimport find_root, join_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. + + """ + cdef np.int_t root = n + while (forest[root] < root): + root = forest[root] + return root + +cdef set_root(np.int_t *forest, np.int_t n, np.int_t root): + """ + Set all nodes on a path to point to new_root. + + """ + cdef np.int_t j + while (forest[n] < n): + j = forest[n] + forest[n] = root + n = j + + forest[n] = root + + +cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): + """Join two trees containing nodes n and m. + + """ + cdef np.int_t root = find_root(forest, n) + cdef np.int_t root_m + + if (n != m): + root_m = find_root(forest, m) + + if (root > root_m): + root = root_m + + set_root(forest, n, root) + set_root(forest, m, root) def felzenszwalb_segmentation(image, k, sigma=0.8): @@ -37,7 +79,7 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): edge_queue = np.argsort(costs) cdef np.int_t *segments_p = segments.data cdef np.int_t seg_new - segment_size = defaultdict(lambda: 1) + cdef np.ndarray[np.int_t, ndim=1] segment_size = np.ones(width * height, dtype=np.int) # inner cost of segments cint = defaultdict(lambda: 0) for edge, cost in zip(edges[edge_queue], costs[edge_queue]): @@ -53,4 +95,8 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): seg_new = find_root(segments_p, seg0) segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] cint[seg_new] = cost - return segments + # unravel the union find tree + old = np.zeros_like(flat) + while (old != flat).any(): + old = flat + flat = flat[flat] From 888d176034369ad47b7333366a6a559228bdef8c Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 20:56:50 +0200 Subject: [PATCH 14/63] ENH much faster. --- skimage/segmentation/felzenszwalb.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/felzenszwalb.pyx index 11f3e93b..8a927fa0 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/felzenszwalb.pyx @@ -78,10 +78,11 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): # and inner cost, then start greedy iteration over edges. edge_queue = np.argsort(costs) cdef np.int_t *segments_p = segments.data - cdef np.int_t seg_new cdef np.ndarray[np.int_t, ndim=1] segment_size = np.ones(width * height, dtype=np.int) # inner cost of segments - cint = defaultdict(lambda: 0) + cdef np.ndarray[np.float_t, ndim=1] cint = np.zeros(width * height) + cdef int seg0, seg1, seg_new + cdef float cost, inner_cost0, inner_cost1 for edge, cost in zip(edges[edge_queue], costs[edge_queue]): seg0 = find_root(segments_p, edge[0]) seg1 = find_root(segments_p, edge[1]) @@ -100,3 +101,4 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): while (old != flat).any(): old = flat flat = flat[flat] + return flat.reshape((width, height)) From 461d4be549b051dda2ceac8582ec69663c636740 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 20:57:10 +0200 Subject: [PATCH 15/63] forgot a line :-/ --- skimage/segmentation/felzenszwalb.pyx | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/felzenszwalb.pyx index 8a927fa0..6b3fd90a 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/felzenszwalb.pyx @@ -97,6 +97,7 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] cint[seg_new] = cost # unravel the union find tree + flat = segments.ravel() old = np.zeros_like(flat) while (old != flat).any(): old = flat From 0c198998256138be9bb034b1203d320bb3734c9e Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 21:26:45 +0200 Subject: [PATCH 16/63] enh cythonizing some arrays --- skimage/segmentation/felzenszwalb.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/felzenszwalb.pyx index 6b3fd90a..1275fa6a 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/felzenszwalb.pyx @@ -64,8 +64,8 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): 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:])) - costs = np.hstack([right_cost.ravel(), down_cost.ravel(), - dright_cost.ravel(), uright_cost.ravel()]) + 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: width, height = image.shape[:2] cdef np.ndarray[np.int_t, ndim=2] segments = np.arange(width * height).reshape(width, height) @@ -73,7 +73,7 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): 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()] - edges = np.vstack([right_edges, down_edges, dright_edges, uright_edges]) + 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) From 7a5e7e49eac35894a86c05099d7af2dc27d9b616 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 21:39:31 +0200 Subject: [PATCH 17/63] ENH reasonable speed for felzenszwalbs's segmentation --- skimage/segmentation/felzenszwalb.pyx | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/felzenszwalb.pyx index 1275fa6a..d45adeb3 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/felzenszwalb.pyx @@ -3,6 +3,7 @@ cimport numpy as np from collections import defaultdict import scipy + #from ..util import img_as_float #from ..color import rgb2grey #from skimage.morphology.ccomp cimport find_root, join_trees @@ -58,8 +59,6 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): image = scipy.ndimage.gaussian_filter(image, sigma=sigma) # compute edge weights in 8 connectivity: - #right_cost = np.sum((image[1:, :, :] - image[:-1, :, :]) ** 2, axis=2) - #down_cost = np.sum((image[:, 1:, :] - image[:, :-1, :]) ** 2, axis=2) 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])) @@ -77,25 +76,35 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): # 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 = segments.data + cdef np.int_t *edges_p = edges.data + cdef np.float_t *costs_p = 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 cdef float cost, inner_cost0, inner_cost1 - for edge, cost in zip(edges[edge_queue], costs[edge_queue]): - seg0 = find_root(segments_p, edge[0]) - seg1 = find_root(segments_p, edge[1]) + # set costs_p back one. we increase it before we use it + # since we might continue before that. + costs_p -= 1 + for e in xrange(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] + k / segment_size[seg0] inner_cost1 = cint[seg1] + k / segment_size[seg1] - if cost < min(inner_cost0, inner_cost1): + 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] = cost + cint[seg_new] = costs_p[0] + # unravel the union find tree flat = segments.ravel() old = np.zeros_like(flat) From 07fb8d0c03cd9586feeeebddf2b3f5450d4254e1 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 22:23:28 +0200 Subject: [PATCH 18/63] ENH felzenszwalb for color images --- .../plot_felzenszwalb_segmentation.py | 24 +++++-- .../{felzenszwalb.pyx => _felzenszwalb.pyx} | 33 +++++++-- skimage/segmentation/felzenszwalb.py | 67 +++++++++++++++++++ skimage/segmentation/setup.py | 4 +- 4 files changed, 114 insertions(+), 14 deletions(-) rename skimage/segmentation/{felzenszwalb.pyx => _felzenszwalb.pyx} (81%) create mode 100644 skimage/segmentation/felzenszwalb.py diff --git a/doc/examples/plot_felzenszwalb_segmentation.py b/doc/examples/plot_felzenszwalb_segmentation.py index dbcb2abb..0b8a2e8c 100644 --- a/doc/examples/plot_felzenszwalb_segmentation.py +++ b/doc/examples/plot_felzenszwalb_segmentation.py @@ -3,11 +3,25 @@ import numpy as np from skimage.data import lena from skimage.segmentation import felzenszwalb_segmentation +from skimage.util import img_as_float -img = lena() -segments = felzenszwalb_segmentation(img, k=1000) -plt.imshow(img) -plt.figure() -plt.imshow(segments) +img = img_as_float(lena()) +segments = felzenszwalb_segmentation(img, scale=1) +segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) + +plt.subplot(131, title="original") +plt.imshow(img, interpolation='nearest') + +plt.subplot(132, title="superpixels") +# shuffle the labels for better visualization +permuted_labels = np.random.permutation(segments.max() + 1) +plt.imshow(permuted_labels[segments], interpolation='nearest') + +plt.subplot(133, title="mean color") +colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in + xrange(img.shape[2])] +counts = np.bincount(segments.ravel()) +colors = np.vstack(colors) / counts +plt.imshow(colors.T[segments], interpolation='nearest') plt.show() print("num segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx similarity index 81% rename from skimage/segmentation/felzenszwalb.pyx rename to skimage/segmentation/_felzenszwalb.pyx index d45adeb3..cbdaac89 100644 --- a/skimage/segmentation/felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -51,11 +51,30 @@ cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): set_root(forest, m, root) -def felzenszwalb_segmentation(image, k, sigma=0.8): - k = float(k) - #image = img_as_float(image) - #image = rgb2grey(image) - image = image[:, :, 0] +def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): + """Computes Felsenszwalb's efficient graph based segmentation for a single channel. + + Parameters + ---------- + image: ndarray, [width, height] + Input image + + scale: float + Free parameter. Higher means larger clusters. + For 0-255 data, hundereds are good. + + sigma: float + Width of Gaussian kernel used in preprocessing. + + Returns + ------- + segment_mask: ndarray, [width, height] + 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)) + scale = float(scale) image = scipy.ndimage.gaussian_filter(image, sigma=sigma) # compute edge weights in 8 connectivity: @@ -96,8 +115,8 @@ def felzenszwalb_segmentation(image, k, sigma=0.8): costs_p += 1 if seg0 == seg1: continue - inner_cost0 = cint[seg0] + k / segment_size[seg0] - inner_cost1 = cint[seg1] + k / segment_size[seg1] + 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) diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py new file mode 100644 index 00000000..6f2f89fc --- /dev/null +++ b/skimage/segmentation/felzenszwalb.py @@ -0,0 +1,67 @@ +import warnings +import numpy as np + +from ._felzenszwalb import felzenszwalb_segmentation_grey + +from IPython.core.debugger import Tracer +tracer = Tracer() + + +def felzenszwalb_segmentation(image, scale=200, sigma=0.8): + """Computes Felsenszwalb's segmentation for multi channel images. + + Calls the algorithm on each channel separately, then combines + using "and", i.e. two pixels are in the same segment if they are + in the same segment for each channel. + + Parameters + ---------- + image: ndarray, [width, height] + Input image + + scale: float + Free parameter. Higher means larger clusters. + For 0-255 data, hundereds are good. + + sigma: float + Width of Gaussian kernel used in preprocessing. + + Returns + ------- + segment_mask: ndarray, [width, height] + Integer mask indicating segment labels. + """ + + #image = img_as_float(image) + if image.ndim == 2: + # assume single channel image + return felzenszwalb_segmentation_grey(image, scale=scale, sigma=sigma) + + elif image.ndim != 3: + raise ValueError("Got image with ndim=%d, don't know" + " what to do." % 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 xrange(n_channels): + channel = np.ascontiguousarray(image[:, :, c]) + seg = felzenszwalb_segmentation_grey(channel, scale=scale, sigma=sigma) + segmentations.append(seg) + + # put pixels in same segment only if in the same segment in all images + # we do this by combining the channels to one number + segmentations = [np.unique(s, return_inverse=True)[1] for s in + segmentations] + n0 = max(segmentations[0]) + n1 = max(segmentations[1]) + hasher = np.array([n1 * n0, n0, 1]) + segmentations = np.dstack(segmentations).reshape(-1, n_channels) + segmentation = np.dot(segmentations, hasher) + # make segment labels consecutive numbers starting at 0 + labels = np.unique(segmentation, return_inverse=True)[1] + return labels.reshape(image.shape[:2]) diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index a197555f..18713881 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -11,8 +11,8 @@ def configuration(parent_package='', top_path=None): config = Configuration('segmentation', parent_package, top_path) - cython(['felzenszwalb.pyx'], working_path=base_path) - config.add_extension('felzenszwalb', sources=['felzenszwalb.c'], + cython(['_felzenszwalb.pyx'], working_path=base_path) + config.add_extension('_felzenszwalb', sources=['_felzenszwalb.c'], include_dirs=[get_numpy_include_dirs()]) cython(['quickshift.pyx'], working_path=base_path) config.add_extension('quickshift', sources=['quickshift.c'], From 80b439bb4aa2fbdd8c7590f064acc43805919a3d Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 22:47:29 +0200 Subject: [PATCH 19/63] ENH Polish examples. --- .../plot_felzenszwalb_segmentation.py | 20 ++++++++++++- doc/examples/plot_quickshift.py | 30 +++++++++++++++---- 2 files changed, 44 insertions(+), 6 deletions(-) diff --git a/doc/examples/plot_felzenszwalb_segmentation.py b/doc/examples/plot_felzenszwalb_segmentation.py index 0b8a2e8c..27581fb1 100644 --- a/doc/examples/plot_felzenszwalb_segmentation.py +++ b/doc/examples/plot_felzenszwalb_segmentation.py @@ -1,3 +1,21 @@ +""" +================================================= +Felzenszwalb's efficient graph based segmentation +================================================= + +This fast 2d image segmentation algorithm, proposed in [1]_ is popular in the +computer vision community. It is often used to extract "superpixels", small +homogeneous image regions, which build the basis for further processing. + +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 +""" +print __doc__ + import matplotlib.pyplot as plt import numpy as np @@ -23,5 +41,5 @@ colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in counts = np.bincount(segments.ravel()) colors = np.vstack(colors) / counts plt.imshow(colors.T[segments], interpolation='nearest') +print("number of segments: %d" % len(np.unique(segments))) plt.show() -print("num segments: %d" % len(np.unique(segments))) diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index b21c56e8..29f44899 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -1,3 +1,26 @@ +""" +============================= +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 color+coordinate space, +see [1]_ It is often used to extract "superpixels", small homogeneous image +regions, which build the basis for further processing. + +One of the benefits of quickshift is that it actually computes a +hierarchical segmentation on multiple scales simultaneously. + +Quickshift has two parameters, one controlling the scale of the local +density approximation, the other selecting a level in the hierarchical +segmentation that is produced. + +.. [1] Quick shift and kernel methods for mode seeking, Vedaldi, A. and Soatto, S. + European Conference on Computer Vision, 2008 +""" +print __doc__ + import matplotlib.pyplot as plt import numpy as np @@ -5,11 +28,8 @@ from skimage.data import lena from skimage.segmentation import quickshift from skimage.util import img_as_float -from IPython.core.debugger import Tracer -tracer = Tracer() - -img = img_as_float(lena()) +img = img_as_float(lena())[::2, ::2, :].copy("C") segments = quickshift(img, sigma=5, tau=20) segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) @@ -27,5 +47,5 @@ colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in counts = np.bincount(segments.ravel()) colors = np.vstack(colors) / counts plt.imshow(colors.T[segments], interpolation='nearest') +print("number of segments: %d" % len(np.unique(segments))) plt.show() -print("num segments: %d" % len(np.unique(segments))) From 83616f0254bf33db6f886cdce2ebe1602c356d0a Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 23:01:24 +0200 Subject: [PATCH 20/63] DOC more docs.... --- skimage/segmentation/felzenszwalb.py | 18 +++++++++++++++--- skimage/segmentation/quickshift.pyx | 23 ++++++++++++++++++----- 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py index 6f2f89fc..b0467d57 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.py @@ -3,13 +3,20 @@ import numpy as np from ._felzenszwalb import felzenszwalb_segmentation_grey -from IPython.core.debugger import Tracer -tracer = Tracer() - def felzenszwalb_segmentation(image, scale=200, sigma=0.8): """Computes Felsenszwalb's segmentation for multi channel images. + 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. + Calls the algorithm on each channel separately, then combines using "and", i.e. two pixels are in the same segment if they are in the same segment for each channel. @@ -30,6 +37,11 @@ def felzenszwalb_segmentation(image, scale=200, sigma=0.8): ------- segment_mask: ndarray, [width, height] 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) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 7e27fc56..b47e5caa 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -10,7 +10,9 @@ cdef extern from "math.h": def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, tau=10, return_tree=False): - """Computes quickshift clustering in RGB-(x,y) space. + """Segments image using quickshift clustering in Color-(x,y) space. + + Produces an oversegmentation of the image using the quickshift mode-seeking algorithm. Parameters ---------- @@ -29,26 +31,37 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta ------- segment_mask: ndarray, [width, height] Integer mask indicating segment labels. + + Notes + ----- + The authors advocate to convert the image to Lab color space prior to segmentation. + + References + ---------- + .. [1] Quick shift and kernel methods for mode seeking, Vedaldi, A. and Soatto, S. + European Conference on Computer Vision, 2008 + + """ # We compute the distances twice since otherwise - # we might get crazy memory overhead (width * height * windowsize**2) - # if you want to speed up things: computing exp in C is the bottleneck ;) + # we get crazy memory overhead (width * height * windowsize**2) # TODO do smoothing beforehand? # TODO manage borders somehow? + # TODO join orphant roots? # window size for neighboring pixels to consider if sigma < 1: raise ValueError("Sigma should be >= 1") cdef int w = int(2 * sigma) - + cdef int width = image.shape[0] cdef int height = image.shape[1] cdef int channels = image.shape[2] cdef float closest, dist cdef int x, y, xx, yy, x_, y_ - + cdef np.float_t* image_p = image.data cdef np.float_t* current_pixel_p = image_p cdef np.float_t* current_entry_p From 9a8cb483c4ee649fe0b8ef3f1ae684ec6f05f6aa Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 23:06:45 +0200 Subject: [PATCH 21/63] misc remove profiling outputs from quickshift --- skimage/segmentation/quickshift.pyx | 8 -------- 1 file changed, 8 deletions(-) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index b47e5caa..df4c14b3 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -3,7 +3,6 @@ cimport numpy as np from itertools import product -from time import time cdef extern from "math.h": double exp(double) @@ -67,7 +66,6 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta cdef np.float_t* current_entry_p cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) - start = time() # compute densities for x, y in product(xrange(width), xrange(height)): for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): @@ -81,15 +79,12 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta densities[x, y] += exp(-dist / sigma) current_pixel_p += channels - print("densities: %f" % (time() - start)) - # this will break ties that otherwise would give us headache densities += np.random.normal(scale=0.00001, size=(width, height)) # default parent to self: cdef np.ndarray[dtype=np.int_t, ndim=2] parent = np.arange(width * height).reshape(width, height) cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) - start = time() # find nearest node with higher density current_pixel_p = image_p for x, y in product(xrange(width), xrange(height)): @@ -108,9 +103,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta parent[x, y] = x_ * width + y_ dist_parent[x, y] = closest current_pixel_p += channels - print("parents: %f" % (time() - start)) - start = time() dist_parent_flat = dist_parent.ravel() flat = parent.ravel() flat[dist_parent_flat > tau] = np.arange(width * height)[dist_parent_flat > tau] @@ -118,7 +111,6 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta while (old != flat).any(): old = flat flat = flat[flat] - print("rest: %f" % (time() - start)) flat = flat.reshape(width, height) if return_tree: return flat, parent From 4d10749a0ec979da8a1d086dbc8314486b5d951f Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 17 Jun 2012 23:17:29 +0200 Subject: [PATCH 22/63] DOC document and export felzenszwalb_segmentation_grey, prettify plots for the web. --- doc/examples/plot_felzenszwalb_segmentation.py | 11 +++++++++-- doc/examples/plot_quickshift.py | 9 ++++++++- skimage/segmentation/__init__.py | 4 +++- skimage/segmentation/_felzenszwalb.pyx | 10 ++++++++++ 4 files changed, 30 insertions(+), 4 deletions(-) diff --git a/doc/examples/plot_felzenszwalb_segmentation.py b/doc/examples/plot_felzenszwalb_segmentation.py index 27581fb1..3d06dd2a 100644 --- a/doc/examples/plot_felzenszwalb_segmentation.py +++ b/doc/examples/plot_felzenszwalb_segmentation.py @@ -27,13 +27,17 @@ img = img_as_float(lena()) segments = felzenszwalb_segmentation(img, scale=1) segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) +print("number of segments: %d" % len(np.unique(segments))) + plt.subplot(131, title="original") plt.imshow(img, interpolation='nearest') +plt.axis("off") -plt.subplot(132, title="superpixels") +plt.subplot(132, title="segmentation") # shuffle the labels for better visualization permuted_labels = np.random.permutation(segments.max() + 1) plt.imshow(permuted_labels[segments], interpolation='nearest') +plt.axis("off") plt.subplot(133, title="mean color") colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in @@ -41,5 +45,8 @@ colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in counts = np.bincount(segments.ravel()) colors = np.vstack(colors) / counts plt.imshow(colors.T[segments], interpolation='nearest') -print("number of segments: %d" % len(np.unique(segments))) +plt.axis("off") + +plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, + bottom=0.02, left=0.02, right=0.98) plt.show() diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index 29f44899..d807ae14 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -33,13 +33,17 @@ img = img_as_float(lena())[::2, ::2, :].copy("C") segments = quickshift(img, sigma=5, tau=20) segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) +print("number of segments: %d" % len(np.unique(segments))) + plt.subplot(131, title="original") plt.imshow(img, interpolation='nearest') +plt.axis("off") plt.subplot(132, title="superpixels") # shuffle the labels for better visualization permuted_labels = np.random.permutation(segments.max() + 1) plt.imshow(permuted_labels[segments], interpolation='nearest') +plt.axis("off") plt.subplot(133, title="mean color") colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in @@ -47,5 +51,8 @@ colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in counts = np.bincount(segments.ravel()) colors = np.vstack(colors) / counts plt.imshow(colors.T[segments], interpolation='nearest') -print("number of segments: %d" % len(np.unique(segments))) +plt.axis("off") + +plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, + bottom=0.02, left=0.02, right=0.98) plt.show() diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 8fa8dfb8..36595240 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,5 +1,7 @@ from .random_walker_segmentation import random_walker from .felzenszwalb import felzenszwalb_segmentation +from .felzenszwalb import felzenszwalb_segmentation_grey from .quickshift import quickshift -__all__ = [random_walker, quickshift, felzenszwalb_segmentation] +__all__ = [random_walker, quickshift, felzenszwalb_segmentation, + felzenszwalb_segmentation_grey] diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index cbdaac89..5a6db993 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -54,6 +54,16 @@ cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): """Computes Felsenszwalb'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 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. + Parameters ---------- image: ndarray, [width, height] From ce26467ad4bba4404d8d2edc46eb3efd06500a40 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 18 Jun 2012 00:37:49 +0200 Subject: [PATCH 23/63] ENH: make quickshift more tolerant to input type, just convert to float. Also keep track of random seed for reproducable tests. Finally, do a unique on the output and add testing. --- doc/examples/plot_quickshift.py | 1 - skimage/segmentation/quickshift.pyx | 28 ++++++++--- skimage/segmentation/tests/test_quickshift.py | 50 +++++++++++++++++++ 3 files changed, 70 insertions(+), 9 deletions(-) create mode 100644 skimage/segmentation/tests/test_quickshift.py diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index d807ae14..f0ace2d4 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -31,7 +31,6 @@ from skimage.util import img_as_float img = img_as_float(lena())[::2, ::2, :].copy("C") segments = quickshift(img, sigma=5, tau=20) -segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) print("number of segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index df4c14b3..35293ea6 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -3,12 +3,14 @@ cimport numpy as np from itertools import product +from ..util import img_as_float + cdef extern from "math.h": double exp(double) -def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, tau=10, return_tree=False): +def quickshift(image, sigma=5, tau=10, return_tree=False, 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. @@ -25,6 +27,8 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta Higher means less clusters. return_tree: bool Whether to return the full segmentation hierarchy tree + random_seed: None or int + Random seed used for breaking ties Returns ------- @@ -42,6 +46,13 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta """ + image = np.atleast_3d(image) + cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = img_as_float(np.ascontiguousarray(image)) + + if random_seed is None: + random_state = np.random.RandomState() + else: + random_state = np.random.RandomState(random_seed) # We compute the distances twice since otherwise # we get crazy memory overhead (width * height * windowsize**2) @@ -55,13 +66,13 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta raise ValueError("Sigma should be >= 1") cdef int w = int(2 * sigma) - cdef int width = image.shape[0] - cdef int height = image.shape[1] - cdef int channels = image.shape[2] + cdef int width = image_c.shape[0] + cdef int height = image_c.shape[1] + cdef int channels = image_c.shape[2] cdef float closest, dist cdef int x, y, xx, yy, x_, y_ - cdef np.float_t* image_p = image.data + cdef np.float_t* image_p = image_c.data cdef np.float_t* current_pixel_p = image_p cdef np.float_t* current_entry_p @@ -74,14 +85,14 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta dist = 0 current_entry_p = current_pixel_p for c in xrange(channels): - dist += (current_pixel_p[c] - image[x_, y_, c])**2 + dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 densities[x, y] += exp(-dist / sigma) current_pixel_p += channels # this will break ties that otherwise would give us headache - densities += np.random.normal(scale=0.00001, size=(width, height)) + densities += random_state.normal(scale=0.00001, size=(width, height)) # default parent to self: cdef np.ndarray[dtype=np.int_t, ndim=2] parent = np.arange(width * height).reshape(width, height) cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) @@ -96,7 +107,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta if densities[x_, y_] > current_density: dist = 0 for c in xrange(channels): - dist += (current_pixel_p[c] - image[x_, y_, c])**2 + dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 if dist < closest: closest = dist @@ -111,6 +122,7 @@ def quickshift(np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image, sigma=5, ta while (old != flat).any(): old = flat flat = flat[flat] + flat = np.unique(flat, return_inverse=True)[1] flat = flat.reshape(width, height) if return_tree: return flat, parent diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py new file mode 100644 index 00000000..b1d17837 --- /dev/null +++ b/skimage/segmentation/tests/test_quickshift.py @@ -0,0 +1,50 @@ +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, 20)) + 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, random_seed=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(): + rnd = np.random.RandomState(0) + img = np.zeros((20, 20, 3)) + img[:10, :10, 0] = 1 + img[10:, :10, 1] = 1 + img[10:, 10:, 2] = 1 + img += 0.2 * rnd.normal(size=img.shape) + img[img > 1] = 1 + img[img < 0] = 0 + seg = quickshift(img, random_seed=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, sigma=1, tau=3, random_seed=0) + # very oversegmented: + assert_equal(len(np.unique(seg2)), 30) + # 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() From 8fa5427afc7d888e9080982401b1cd0376299256 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 18 Jun 2012 08:26:12 +0200 Subject: [PATCH 24/63] FIX build problem and cython problem resolved. --- skimage/morphology/ccomp.pyx | 1 - skimage/segmentation/_felzenszwalb.pyx | 47 +------------------------- 2 files changed, 1 insertion(+), 47 deletions(-) diff --git a/skimage/morphology/ccomp.pyx b/skimage/morphology/ccomp.pyx index a1d5f303..ffb2cf11 100644 --- a/skimage/morphology/ccomp.pyx +++ b/skimage/morphology/ccomp.pyx @@ -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. diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index 5a6db993..afa39ad8 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -3,52 +3,7 @@ cimport numpy as np from collections import defaultdict import scipy - -#from ..util import img_as_float -#from ..color import rgb2grey -#from skimage.morphology.ccomp cimport find_root, join_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. - - """ - cdef np.int_t root = n - while (forest[root] < root): - root = forest[root] - return root - -cdef set_root(np.int_t *forest, np.int_t n, np.int_t root): - """ - Set all nodes on a path to point to new_root. - - """ - cdef np.int_t j - while (forest[n] < n): - j = forest[n] - forest[n] = root - n = j - - forest[n] = root - - -cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): - """Join two trees containing nodes n and m. - - """ - cdef np.int_t root = find_root(forest, n) - cdef np.int_t root_m - - if (n != m): - root_m = find_root(forest, m) - - if (root > root_m): - root = root_m - - set_root(forest, n, root) - set_root(forest, m, root) +from skimage.morphology.ccomp cimport find_root, join_trees def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): From 08df2a5103b8ab038ceee654930dddb0cb68fc09 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Tue, 19 Jun 2012 21:22:48 +0200 Subject: [PATCH 25/63] enh: minor simplifications --- skimage/segmentation/quickshift.pyx | 46 +++++++++++++++-------------- 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 35293ea6..2b6bdf24 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -1,5 +1,6 @@ import numpy as np cimport numpy as np +cimport cython from itertools import product @@ -10,6 +11,9 @@ cdef extern from "math.h": double exp(double) +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): """Segments image using quickshift clustering in Color-(x,y) space. @@ -70,24 +74,22 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): cdef int height = image_c.shape[1] cdef int channels = image_c.shape[2] cdef float closest, dist - cdef int x, y, xx, yy, x_, y_ + cdef int x, y, x_, y_ cdef np.float_t* image_p = image_c.data cdef np.float_t* current_pixel_p = image_p - cdef np.float_t* current_entry_p cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) # compute densities for x, y in product(xrange(width), xrange(height)): - for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): - x_, y_ = x + xx, y + yy - if 0 <= x_ < width and 0 <= y_ < height: - dist = 0 - current_entry_p = current_pixel_p - for c in xrange(channels): - dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 - dist += (x - x_)**2 + (y - y_)**2 - densities[x, y] += exp(-dist / sigma) + x_min, x_max = max(x - w, 0), min(x + w + 1, width) + y_min, y_max = max(y - w, 0), min(y + w + 1, height) + for x_, y_ in product(xrange(x_min, x_max), xrange(y_min, y_max)): + dist = 0 + for c in xrange(channels): + dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 + dist += (x - x_)**2 + (y - y_)**2 + densities[x, y] += exp(-dist / sigma) current_pixel_p += channels # this will break ties that otherwise would give us headache @@ -101,17 +103,17 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): for x, y in product(xrange(width), xrange(height)): current_density = densities[x, y] closest = np.inf - for xx, yy in product(xrange(-w / 2, w / 2 + 1), repeat=2): - x_, y_ = x + xx, y + yy - if 0 <= x_ < width and 0 <= y_ < height: - if densities[x_, y_] > current_density: - dist = 0 - for c in xrange(channels): - dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 - dist += (x - x_)**2 + (y - y_)**2 - if dist < closest: - closest = dist - parent[x, y] = x_ * width + y_ + x_min, x_max = max(x - w, 0), min(x + w + 1, width) + y_min, y_max = max(y - w, 0), min(y + w + 1, height) + for x_, y_ in product(xrange(x_min, x_max), xrange(y_min, y_max)): + if densities[x_, y_] > current_density: + dist = 0 + for c in xrange(channels): + dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 + dist += (x - x_)**2 + (y - y_)**2 + if dist < closest: + closest = dist + parent[x, y] = x_ * width + y_ dist_parent[x, y] = closest current_pixel_p += channels From f0a7212c4f710fb3798ecb3c1b8e6fe8e910495c Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Tue, 19 Jun 2012 21:56:51 +0200 Subject: [PATCH 26/63] ENH Rename parameters in quickshift, add "ratio" --- doc/examples/plot_quickshift.py | 3 +-- skimage/segmentation/quickshift.pyx | 19 +++++++++++-------- skimage/segmentation/tests/test_quickshift.py | 4 ++-- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index f0ace2d4..1285ad28 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -28,9 +28,8 @@ from skimage.data import lena from skimage.segmentation import quickshift from skimage.util import img_as_float - img = img_as_float(lena())[::2, ::2, :].copy("C") -segments = quickshift(img, sigma=5, tau=20) +segments = quickshift(img, kernel_size=5, max_dist=20) print("number of segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 2b6bdf24..1822055b 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -14,7 +14,7 @@ cdef extern from "math.h": @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): +def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, 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. @@ -23,10 +23,13 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): ---------- image: ndarray, [width, height, channels] Input image - sigma: float + 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 less clusters. - tau: float + max_dist: float Cut-off point for data distances. Higher means less clusters. return_tree: bool @@ -51,7 +54,7 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): """ image = np.atleast_3d(image) - cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = img_as_float(np.ascontiguousarray(image)) + cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = img_as_float(np.ascontiguousarray(image)) * ratio if random_seed is None: random_state = np.random.RandomState() @@ -66,9 +69,9 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): # TODO join orphant roots? # window size for neighboring pixels to consider - if sigma < 1: + if kernel_size < 1: raise ValueError("Sigma should be >= 1") - cdef int w = int(2 * sigma) + cdef int w = int(2 * kernel_size) cdef int width = image_c.shape[0] cdef int height = image_c.shape[1] @@ -89,7 +92,7 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): for c in xrange(channels): dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 - densities[x, y] += exp(-dist / sigma) + densities[x, y] += exp(-dist / kernel_size) current_pixel_p += channels # this will break ties that otherwise would give us headache @@ -119,7 +122,7 @@ def quickshift(image, sigma=5, tau=10, return_tree=False, random_seed=None): dist_parent_flat = dist_parent.ravel() flat = parent.ravel() - flat[dist_parent_flat > tau] = np.arange(width * height)[dist_parent_flat > tau] + flat[dist_parent_flat > max_dist] = np.arange(width * height)[dist_parent_flat > max_dist] old = np.zeros_like(flat) while (old != flat).any(): old = flat diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py index b1d17837..a904837c 100644 --- a/skimage/segmentation/tests/test_quickshift.py +++ b/skimage/segmentation/tests/test_quickshift.py @@ -37,9 +37,9 @@ def test_color(): assert_array_equal(seg[:10, 10:], 1) assert_array_equal(seg[10:, 10:], 2) - seg2 = quickshift(img, sigma=1, tau=3, random_seed=0) + seg2 = quickshift(img, kernel_size=1, max_dist=3, random_seed=0) # very oversegmented: - assert_equal(len(np.unique(seg2)), 30) + assert_equal(len(np.unique(seg2)), 18) # still don't cross lines assert_true((seg2[9, :] != seg2[10, :]).all()) assert_true((seg2[:, 9] != seg2[:, 10]).all()) From a7c98cb67a4eaca083f86a6c6fd1964a90b05b2f Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Tue, 19 Jun 2012 22:32:44 +0200 Subject: [PATCH 27/63] Remove felzenszwalb_segmentation_gray again since it just complicates the interface. --- skimage/segmentation/__init__.py | 4 +--- skimage/segmentation/_felzenszwalb.pyx | 10 ++++++---- skimage/segmentation/felzenszwalb.py | 14 ++++++-------- 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 36595240..8fa8dfb8 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,7 +1,5 @@ from .random_walker_segmentation import random_walker from .felzenszwalb import felzenszwalb_segmentation -from .felzenszwalb import felzenszwalb_segmentation_grey from .quickshift import quickshift -__all__ = [random_walker, quickshift, felzenszwalb_segmentation, - felzenszwalb_segmentation_grey] +__all__ = [random_walker, quickshift, felzenszwalb_segmentation] diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index afa39ad8..d058975d 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -1,12 +1,13 @@ import numpy as np cimport numpy as np -from collections import defaultdict import scipy from skimage.morphology.ccomp cimport find_root, join_trees +from ..util import img_as_float -def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): + +def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): """Computes Felsenszwalb's efficient graph based segmentation for a single channel. Produces an oversegmentation of a 2d image using a fast, minimum spanning @@ -26,7 +27,6 @@ def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): scale: float Free parameter. Higher means larger clusters. - For 0-255 data, hundereds are good. sigma: float Width of Gaussian kernel used in preprocessing. @@ -39,6 +39,7 @@ def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): 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) scale = float(scale) image = scipy.ndimage.gaussian_filter(image, sigma=sigma) @@ -88,11 +89,12 @@ def felzenszwalb_segmentation_grey(image, scale=200, sigma=0.8): seg_new = find_root(segments_p, seg0) segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] cint[seg_new] = costs_p[0] - + # 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((width, height)) diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py index b0467d57..be879e13 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.py @@ -1,11 +1,11 @@ import warnings import numpy as np -from ._felzenszwalb import felzenszwalb_segmentation_grey +from ._felzenszwalb import _felzenszwalb_segmentation_grey -def felzenszwalb_segmentation(image, scale=200, sigma=0.8): - """Computes Felsenszwalb's segmentation for multi channel images. +def felzenszwalb_segmentation(image, scale=1, sigma=0.8): + """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. @@ -47,7 +47,7 @@ def felzenszwalb_segmentation(image, scale=200, sigma=0.8): #image = img_as_float(image) if image.ndim == 2: # assume single channel image - return felzenszwalb_segmentation_grey(image, scale=scale, sigma=sigma) + return _felzenszwalb_segmentation_grey(image, scale=scale, sigma=sigma) elif image.ndim != 3: raise ValueError("Got image with ndim=%d, don't know" @@ -62,13 +62,11 @@ def felzenszwalb_segmentation(image, scale=200, sigma=0.8): # compute quickshift for each channel for c in xrange(n_channels): channel = np.ascontiguousarray(image[:, :, c]) - seg = felzenszwalb_segmentation_grey(channel, scale=scale, sigma=sigma) - segmentations.append(seg) + s = _felzenszwalb_segmentation_grey(channel, scale=scale, sigma=sigma) + 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 - segmentations = [np.unique(s, return_inverse=True)[1] for s in - segmentations] n0 = max(segmentations[0]) n1 = max(segmentations[1]) hasher = np.array([n1 * n0, n0, 1]) From d2e226fe59ae1b9aa43207e10ea84d85ea160e7c Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Tue, 19 Jun 2012 22:33:14 +0200 Subject: [PATCH 28/63] ENH tests for Felzenszwalbs segmentation, fixed off-by-one error --- skimage/segmentation/felzenszwalb.py | 4 +- .../segmentation/tests/test_felzenszwalb.py | 39 +++++++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) create mode 100644 skimage/segmentation/tests/test_felzenszwalb.py diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py index be879e13..54846102 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.py @@ -67,8 +67,8 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8): # 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 = max(segmentations[0]) - n1 = max(segmentations[1]) + n0 = segmentations[0].max() + 1 + n1 = segmentations[1].max() + 1 hasher = np.array([n1 * n0, n0, 1]) segmentations = np.dstack(segmentations).reshape(-1, n_channels) segmentation = np.dot(segmentations, hasher) diff --git a/skimage/segmentation/tests/test_felzenszwalb.py b/skimage/segmentation/tests/test_felzenszwalb.py new file mode 100644 index 00000000..d645ab9d --- /dev/null +++ b/skimage/segmentation/tests/test_felzenszwalb.py @@ -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_segmentation + + +def test_grey(): + # very weak tests. This algorithm is pretty unstable. + img = np.zeros((20, 20)) + img[:10, 10:] = 0.2 + img[10:, :10] = 0.4 + img[10:, 10:] = 0.6 + seg = felzenszwalb_segmentation(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, 20, 3)) + img[:10, :10, 0] = 1 + img[10:, :10, 1] = 1 + img[10:, 10:, 2] = 1 + seg = felzenszwalb_segmentation(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], 3) + assert_array_equal(seg[:10, 10:], 1) + assert_array_equal(seg[10:, 10:], 2) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() From 05cc863f3f2090a1eb10f96ba185cd0c9225d2f8 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Tue, 19 Jun 2012 23:58:40 +0200 Subject: [PATCH 29/63] First draft for numpy based km_segmentation --- skimage/segmentation/__init__.py | 4 ++- skimage/segmentation/km_segmentation.py | 42 +++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 1 deletion(-) create mode 100644 skimage/segmentation/km_segmentation.py diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 8fa8dfb8..2de27ad1 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,5 +1,7 @@ from .random_walker_segmentation import random_walker from .felzenszwalb import felzenszwalb_segmentation +from .km_segmentation import km_segmentation from .quickshift import quickshift -__all__ = [random_walker, quickshift, felzenszwalb_segmentation] +__all__ = [random_walker, quickshift, felzenszwalb_segmentation, + km_segmentation] diff --git a/skimage/segmentation/km_segmentation.py b/skimage/segmentation/km_segmentation.py new file mode 100644 index 00000000..fa307ab3 --- /dev/null +++ b/skimage/segmentation/km_segmentation.py @@ -0,0 +1,42 @@ +import numpy as np + + +def km_segmentation(image, n_segments=100, ratio=50, max_iter=100): + # initialize on grid: + height, width = image.shape[:2] + # approximate grid size for desired n_segments + step = 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 = image[means_y, means_x, :] + means = np.dstack([means_y, means_x, means_color]).reshape(-1, 5) + image = np.dstack([grid_y, grid_x, image * ratio]) + + nearest_mean = np.zeros((height, width), dtype=np.int) + distance = np.ones((height, width), dtype=np.float) * np.inf + for i in xrange(max_iter): + print("iteration %d" % i) + nearest_mean_old = nearest_mean.copy() + # assign pixels to means + for k, mean in enumerate(means): + # compute windows: + y_min = int(max(mean[0] - 2 * step, 0)) + y_max = int(min(mean[0] + 2 * step, height)) + x_min = int(max(mean[1] - 2 * step, 0)) + x_max = int(min(mean[1] + 2 * step, height)) + search_window = image[y_min:y_max + 1, x_min:x_max + 1] + dist_mean = np.sum((search_window - mean) ** 2, axis=2) + assign = distance[y_min:y_max + 1, x_min:x_max + 1] > dist_mean + nearest_mean[y_min:y_max + 1, x_min:x_max + 1][assign] = k + distance[y_min:y_max + 1, x_min:x_max + 1][assign] = \ + dist_mean[assign] + if (nearest_mean == nearest_mean_old).all(): + break + # recompute means: + means = [np.bincount(nearest_mean.ravel(), image[:, :, j].ravel()) + for j in xrange(5)] + in_mean = np.bincount(nearest_mean.ravel()) + means = (np.vstack(means) / in_mean).T + return nearest_mean From ccfb89b9570598f41cff7ad4d34f3236c1d911f3 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 22 Jun 2012 23:08:59 +0200 Subject: [PATCH 30/63] FIX Tried to address @stefanv's comments on the PR. --- .../plot_felzenszwalb_segmentation.py | 25 +++++++++---------- doc/examples/plot_quickshift.py | 24 ++++++++---------- skimage/segmentation/__init__.py | 2 +- skimage/segmentation/_felzenszwalb.pyx | 6 ++--- skimage/segmentation/felzenszwalb.py | 10 +++----- skimage/segmentation/quickshift.pyx | 4 +-- 6 files changed, 31 insertions(+), 40 deletions(-) diff --git a/doc/examples/plot_felzenszwalb_segmentation.py b/doc/examples/plot_felzenszwalb_segmentation.py index 3d06dd2a..18ac0f80 100644 --- a/doc/examples/plot_felzenszwalb_segmentation.py +++ b/doc/examples/plot_felzenszwalb_segmentation.py @@ -29,24 +29,23 @@ segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) print("number of segments: %d" % len(np.unique(segments))) -plt.subplot(131, title="original") -plt.imshow(img, interpolation='nearest') -plt.axis("off") -plt.subplot(132, title="segmentation") -# shuffle the labels for better visualization -permuted_labels = np.random.permutation(segments.max() + 1) -plt.imshow(permuted_labels[segments], interpolation='nearest') -plt.axis("off") +fig, (ax_org, ax_sp, ax_mean) = plt.subplots(1, 3) +ax_org.set_title("original") +ax_org.imshow(img, interpolation='nearest') +ax_org.axis("off") + +ax_sp.set_title("superpixels") +ax_sp.imshow(segments, interpolation='nearest', cmap=plt.cm.prism) +ax_sp.axis("off") -plt.subplot(133, title="mean color") colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in xrange(img.shape[2])] counts = np.bincount(segments.ravel()) colors = np.vstack(colors) / counts -plt.imshow(colors.T[segments], interpolation='nearest') -plt.axis("off") - -plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, +ax_mean.set_title("mean color") +ax_mean.imshow(colors.T[segments], interpolation='nearest') +ax_mean.axis("off") +fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0.02, left=0.02, right=0.98) plt.show() diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index 1285ad28..5adf1969 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -33,24 +33,22 @@ segments = quickshift(img, kernel_size=5, max_dist=20) print("number of segments: %d" % len(np.unique(segments))) -plt.subplot(131, title="original") -plt.imshow(img, interpolation='nearest') -plt.axis("off") +fig, (ax_org, ax_sp, ax_mean) = plt.subplots(1, 3) +ax_org.set_title("original") +ax_org.imshow(img, interpolation='nearest') +ax_org.axis("off") -plt.subplot(132, title="superpixels") -# shuffle the labels for better visualization -permuted_labels = np.random.permutation(segments.max() + 1) -plt.imshow(permuted_labels[segments], interpolation='nearest') -plt.axis("off") +ax_sp.set_title("superpixels") +ax_sp.imshow(segments, interpolation='nearest', cmap=plt.cm.prism) +ax_sp.axis("off") -plt.subplot(133, title="mean color") colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in xrange(img.shape[2])] counts = np.bincount(segments.ravel()) colors = np.vstack(colors) / counts -plt.imshow(colors.T[segments], interpolation='nearest') -plt.axis("off") - -plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, +ax_mean.set_title("mean color") +ax_mean.imshow(colors.T[segments], interpolation='nearest') +ax_mean.axis("off") +fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0.02, left=0.02, right=0.98) plt.show() diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 2de27ad1..b2c97448 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -4,4 +4,4 @@ from .km_segmentation import km_segmentation from .quickshift import quickshift __all__ = [random_walker, quickshift, felzenszwalb_segmentation, - km_segmentation] + km_segmentation] diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index d058975d..a677020a 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -22,12 +22,10 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): Parameters ---------- - image: ndarray, [width, height] + image: (width, height) ndarray Input image - scale: float Free parameter. Higher means larger clusters. - sigma: float Width of Gaussian kernel used in preprocessing. @@ -74,7 +72,7 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): # set costs_p back one. we increase it before we use it # since we might continue before that. costs_p -= 1 - for e in xrange(costs.size): + 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 diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py index 54846102..cfcf4f23 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.py @@ -23,13 +23,10 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8): Parameters ---------- - image: ndarray, [width, height] + image: (width, height) ndarray Input image - scale: float Free parameter. Higher means larger clusters. - For 0-255 data, hundereds are good. - sigma: float Width of Gaussian kernel used in preprocessing. @@ -69,9 +66,8 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8): # we do this by combining the channels to one number n0 = segmentations[0].max() + 1 n1 = segmentations[1].max() + 1 - hasher = np.array([n1 * n0, n0, 1]) - segmentations = np.dstack(segmentations).reshape(-1, n_channels) - segmentation = np.dot(segmentations, hasher) + 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]) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 1822055b..b0a9b2d6 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -21,7 +21,7 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r Parameters ---------- - image: ndarray, [width, height, channels] + image: (width, height, channels) ndarray Input image ratio: float, between 0 and 1. Balances color-space proximity and image-space proximity. @@ -54,7 +54,7 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r """ image = np.atleast_3d(image) - cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = img_as_float(np.ascontiguousarray(image)) * ratio + cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = np.ascontiguousarray(img_as_float(image)) * ratio if random_seed is None: random_state = np.random.RandomState() From 49bc44c6e9686a1f8d3d9a5a62b458c11442264a Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 22 Jun 2012 23:10:58 +0200 Subject: [PATCH 31/63] FIXed test to work with the fixed "hashing" of colors --- skimage/segmentation/tests/test_felzenszwalb.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/tests/test_felzenszwalb.py b/skimage/segmentation/tests/test_felzenszwalb.py index d645ab9d..f6cca31b 100644 --- a/skimage/segmentation/tests/test_felzenszwalb.py +++ b/skimage/segmentation/tests/test_felzenszwalb.py @@ -29,9 +29,9 @@ def test_color(): # 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], 2) assert_array_equal(seg[:10, 10:], 1) - assert_array_equal(seg[10:, 10:], 2) + assert_array_equal(seg[10:, 10:], 3) if __name__ == '__main__': From a779952619d65307a8ec6e29f1e416a2eeab9688 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 22 Jun 2012 23:16:25 +0200 Subject: [PATCH 32/63] DOC updates CONTRIBUTORS.txt --- CONTRIBUTORS.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index baaf064b..b00593c7 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -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. From 7b646ad7eabf4cfe356b327b697643ab62fd1da5 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 23 Jun 2012 01:07:32 +0200 Subject: [PATCH 33/63] fix trying to make this implementation more like slic --- skimage/segmentation/km_segmentation.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/km_segmentation.py b/skimage/segmentation/km_segmentation.py index fa307ab3..1c50a8da 100644 --- a/skimage/segmentation/km_segmentation.py +++ b/skimage/segmentation/km_segmentation.py @@ -1,7 +1,10 @@ import numpy as np +from scipy import ndimage +from ..util import img_as_float -def km_segmentation(image, n_segments=100, ratio=50, max_iter=100): +def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): + image = ndimage.gaussian_filter(img_as_float(image), sigma) # initialize on grid: height, width = image.shape[:2] # approximate grid size for desired n_segments @@ -12,7 +15,11 @@ def km_segmentation(image, n_segments=100, ratio=50, max_iter=100): means_color = image[means_y, means_x, :] means = np.dstack([means_y, means_x, means_color]).reshape(-1, 5) - image = np.dstack([grid_y, grid_x, image * ratio]) + # 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 + print(ratio) + image = np.dstack([grid_y, grid_x, image / ratio]) nearest_mean = np.zeros((height, width), dtype=np.int) distance = np.ones((height, width), dtype=np.float) * np.inf From 8f243667900e48da50069db4ce777c0aa8fb192e Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 2 Jul 2012 21:29:01 +0100 Subject: [PATCH 34/63] starting cython implementation of km_segmentation --- skimage/segmentation/km_segmentation.py | 49 --------------------- skimage/segmentation/km_segmentation.pyx | 55 ++++++++++++++++++++++++ skimage/segmentation/setup.py | 3 ++ 3 files changed, 58 insertions(+), 49 deletions(-) delete mode 100644 skimage/segmentation/km_segmentation.py create mode 100644 skimage/segmentation/km_segmentation.pyx diff --git a/skimage/segmentation/km_segmentation.py b/skimage/segmentation/km_segmentation.py deleted file mode 100644 index 1c50a8da..00000000 --- a/skimage/segmentation/km_segmentation.py +++ /dev/null @@ -1,49 +0,0 @@ -import numpy as np -from scipy import ndimage -from ..util import img_as_float - - -def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): - image = ndimage.gaussian_filter(img_as_float(image), sigma) - # initialize on grid: - height, width = image.shape[:2] - # approximate grid size for desired n_segments - step = 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 = image[means_y, means_x, :] - means = np.dstack([means_y, means_x, means_color]).reshape(-1, 5) - # 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 - print(ratio) - image = np.dstack([grid_y, grid_x, image / ratio]) - - nearest_mean = np.zeros((height, width), dtype=np.int) - distance = np.ones((height, width), dtype=np.float) * np.inf - for i in xrange(max_iter): - print("iteration %d" % i) - nearest_mean_old = nearest_mean.copy() - # assign pixels to means - for k, mean in enumerate(means): - # compute windows: - y_min = int(max(mean[0] - 2 * step, 0)) - y_max = int(min(mean[0] + 2 * step, height)) - x_min = int(max(mean[1] - 2 * step, 0)) - x_max = int(min(mean[1] + 2 * step, height)) - search_window = image[y_min:y_max + 1, x_min:x_max + 1] - dist_mean = np.sum((search_window - mean) ** 2, axis=2) - assign = distance[y_min:y_max + 1, x_min:x_max + 1] > dist_mean - nearest_mean[y_min:y_max + 1, x_min:x_max + 1][assign] = k - distance[y_min:y_max + 1, x_min:x_max + 1][assign] = \ - dist_mean[assign] - if (nearest_mean == nearest_mean_old).all(): - break - # recompute means: - means = [np.bincount(nearest_mean.ravel(), image[:, :, j].ravel()) - for j in xrange(5)] - in_mean = np.bincount(nearest_mean.ravel()) - means = (np.vstack(means) / in_mean).T - return nearest_mean diff --git a/skimage/segmentation/km_segmentation.pyx b/skimage/segmentation/km_segmentation.pyx new file mode 100644 index 00000000..53237618 --- /dev/null +++ b/skimage/segmentation/km_segmentation.pyx @@ -0,0 +1,55 @@ +import numpy as np +cimport numpy as np +from scipy import ndimage +from ..util import img_as_float + + +def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): + image = ndimage.gaussian_filter(img_as_float(image), sigma) + # initialize on grid: + height, width = image.shape[:2] + # approximate grid size for desired n_segments + step = 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 = image[means_y, means_x, :] + cdef np.ndarray[dtype=np.float_t, ndim=2] means = np.dstack([means_y, means_x, means_color]).reshape(-1, 5) + 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 + print(ratio) + cdef np.ndarray[dtype=np.float_t, ndim=3] image_yx = np.dstack([grid_y, grid_x, image / ratio]) + cdef int i, k, x, y, x_min, x_max, y_min, y_max + cdef float 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.ones((height, width), dtype=np.float) * np.inf + for i in xrange(max_iter): + print("iteration %d" % i) + nearest_mean_old = nearest_mean.copy() + # assign pixels to means + for k in xrange(n_means): + # compute windows: + y_min = int(max(means[k, 0] - 2 * step, 0)) + y_max = int(min(means[k, 0] + 2 * step, height)) + x_min = int(max(means[k, 1] - 2 * step, 0)) + x_max = int(min(means[k, 1] + 2 * step, height)) + for x in xrange(x_min, x_max): + for y in xrange(y_min, y_max): + dist_mean = 0 + for c in range(5): + dist_mean += (image_yx[y, x, c] - means[k, c]) ** 2 + if distance[y, x] > dist_mean: + nearest_mean[y, x] = k + distance[y, x] = dist_mean + if (nearest_mean == nearest_mean_old).all(): + break + # recompute means: + means_list = [np.bincount(nearest_mean.ravel(), image_yx[:, :, j].ravel()) + for j in xrange(5)] + in_mean = np.bincount(nearest_mean.ravel()) + means = (np.vstack(means_list) / in_mean).T + return nearest_mean diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index 18713881..0be6b748 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -17,6 +17,9 @@ def configuration(parent_package='', top_path=None): cython(['quickshift.pyx'], working_path=base_path) config.add_extension('quickshift', sources=['quickshift.c'], include_dirs=[get_numpy_include_dirs()]) + cython(['km_segmentation.pyx'], working_path=base_path) + config.add_extension('km_segmentation', sources=['km_segmentation.c'], + include_dirs=[get_numpy_include_dirs()]) return config From 501a6db8ad048ae25396ab6d4da50de475f29e52 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 2 Jul 2012 22:54:49 +0100 Subject: [PATCH 35/63] ENH speedup, means and image use pointers --- skimage/segmentation/km_segmentation.pyx | 34 +++++++++++++++++------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/skimage/segmentation/km_segmentation.pyx b/skimage/segmentation/km_segmentation.pyx index 53237618..2f1dae5d 100644 --- a/skimage/segmentation/km_segmentation.pyx +++ b/skimage/segmentation/km_segmentation.pyx @@ -16,40 +16,56 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): means_color = image[means_y, means_x, :] 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 print(ratio) - cdef np.ndarray[dtype=np.float_t, ndim=3] image_yx = np.dstack([grid_y, grid_x, image / ratio]) + 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 cdef float 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.ones((height, width), dtype=np.float) * np.inf + cdef np.float_t* image_p = image_yx.data + cdef np.float_t* distance_p = distance.data + cdef np.float_t* current_pixel + cdef float tmp for i in xrange(max_iter): print("iteration %d" % i) nearest_mean_old = nearest_mean.copy() + # we construct a new means every iteration, adjust pointer + current_mean = means.data # assign pixels to means for k in xrange(n_means): # compute windows: - y_min = int(max(means[k, 0] - 2 * step, 0)) - y_max = int(min(means[k, 0] + 2 * step, height)) - x_min = int(max(means[k, 1] - 2 * step, 0)) - x_max = int(min(means[k, 1] + 2 * step, height)) - for x in xrange(x_min, x_max): - for y in xrange(y_min, y_max): + 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, height)) + for y in xrange(y_min, y_max): + current_pixel = &image_p[5 * (y * width + x_min)] + for x in xrange(x_min, x_max): + mean_entry = current_mean dist_mean = 0 for c in range(5): - dist_mean += (image_yx[y, x, c] - means[k, c]) ** 2 + # you would think the compiler can optimize this itself. + # mine can't (with O2) + tmp = current_pixel[0] - mean_entry[0] + current_pixel += 1 + mean_entry += 1 + dist_mean += tmp * tmp if distance[y, x] > dist_mean: nearest_mean[y, x] = k distance[y, x] = dist_mean + current_mean += 5 if (nearest_mean == nearest_mean_old).all(): break # recompute means: means_list = [np.bincount(nearest_mean.ravel(), image_yx[:, :, j].ravel()) for j in xrange(5)] in_mean = np.bincount(nearest_mean.ravel()) - means = (np.vstack(means_list) / in_mean).T + means = (np.vstack(means_list) / in_mean).T.copy("C") return nearest_mean From 1c69adb81783a32c6da234dd748dc48c7d58f063 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Tue, 3 Jul 2012 22:26:52 +0100 Subject: [PATCH 36/63] MISC some simplifications, minor speedup --- skimage/segmentation/km_segmentation.pyx | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/skimage/segmentation/km_segmentation.pyx b/skimage/segmentation/km_segmentation.pyx index 2f1dae5d..f8cfb4e5 100644 --- a/skimage/segmentation/km_segmentation.pyx +++ b/skimage/segmentation/km_segmentation.pyx @@ -24,19 +24,19 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): ratio = (ratio / float(step)) ** 2 print(ratio) 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 - cdef float dist_mean + 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.ones((height, width), dtype=np.float) * np.inf cdef np.float_t* image_p = image_yx.data cdef np.float_t* distance_p = distance.data + cdef np.float_t* current_distance cdef np.float_t* current_pixel - cdef float tmp + cdef double tmp for i in xrange(max_iter): + changes = 0 print("iteration %d" % i) - nearest_mean_old = nearest_mean.copy() - # we construct a new means every iteration, adjust pointer current_mean = means.data # assign pixels to means for k in xrange(n_means): @@ -47,6 +47,7 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): x_max = int(min(current_mean[1] + 2 * step, height)) for y in xrange(y_min, y_max): current_pixel = &image_p[5 * (y * width + x_min)] + current_distance = &distance_p[y * width + x_min] for x in xrange(x_min, x_max): mean_entry = current_mean dist_mean = 0 @@ -54,14 +55,16 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): # you would think the compiler can optimize this itself. # mine can't (with O2) tmp = current_pixel[0] - mean_entry[0] + dist_mean += tmp * tmp current_pixel += 1 mean_entry += 1 - dist_mean += tmp * tmp - if distance[y, x] > dist_mean: + # some precision issue here. Doesnt work if testing ">" + if current_distance[0] - dist_mean > 1e-10: nearest_mean[y, x] = k - distance[y, x] = dist_mean + current_distance[0] = dist_mean + changes += 1 + current_distance += 1 current_mean += 5 - if (nearest_mean == nearest_mean_old).all(): break # recompute means: means_list = [np.bincount(nearest_mean.ravel(), image_yx[:, :, j].ravel()) From d9a22d867bacca4ea297cff286813359edac110c Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Wed, 18 Jul 2012 07:25:39 +0100 Subject: [PATCH 37/63] Documentation, example for km_segmentation --- doc/examples/plot_km_segmentation.py | 36 ++++++++++++++++++++++ skimage/segmentation/km_segmentation.pyx | 39 ++++++++++++++++++++++-- 2 files changed, 73 insertions(+), 2 deletions(-) create mode 100644 doc/examples/plot_km_segmentation.py diff --git a/doc/examples/plot_km_segmentation.py b/doc/examples/plot_km_segmentation.py new file mode 100644 index 00000000..ed8cf983 --- /dev/null +++ b/doc/examples/plot_km_segmentation.py @@ -0,0 +1,36 @@ +""" +""" +print __doc__ + +import matplotlib.pyplot as plt +import numpy as np + +from skimage.data import lena +from skimage.segmentation import km_segmentation +from skimage.util import img_as_float + +img = img_as_float(lena()).copy("C") +segments = km_segmentation(img, ratio=2.0, n_segments=200) + +print("number of segments: %d" % len(np.unique(segments))) + +plt.subplot(131, title="original") +plt.imshow(img, interpolation='nearest') +plt.axis("off") + +plt.subplot(132, title="superpixels") +# shuffle the labels for better visualization +plt.imshow(segments, interpolation='nearest', cmap=plt.cm.prism) +plt.axis("off") + +plt.subplot(133, title="mean color") +colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in + xrange(img.shape[2])] +counts = np.bincount(segments.ravel()) +colors = np.vstack(colors) / counts +plt.imshow(colors.T[segments], interpolation='nearest') +plt.axis("off") + +plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, + bottom=0.02, left=0.02, right=0.98) +plt.show() diff --git a/skimage/segmentation/km_segmentation.pyx b/skimage/segmentation/km_segmentation.pyx index f8cfb4e5..d0674be4 100644 --- a/skimage/segmentation/km_segmentation.pyx +++ b/skimage/segmentation/km_segmentation.pyx @@ -1,11 +1,47 @@ import numpy as np cimport numpy as np +from time import time from scipy import ndimage from ..util import img_as_float def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): + """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. + + Returns + ------- + segment_mask: ndarray, [width, height] + Integer mask indicating segment labels. + + Notes + ----- + The image is smoothed using a Gaussian kernel prior to segmentation. + Best results are achieved if the image is given in Lab color space. + + References + ---------- + .. [1] Slic superpixels, Achanta, R. and Shaji, A. and Smith, K. and Lucchi, + A. and Fua, P. and Suesstrunk, S. + Technical Report 2010 + + """ + 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) + # initialize on grid: height, width = image.shape[:2] # approximate grid size for desired n_segments @@ -22,7 +58,6 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.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 - print(ratio) 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 @@ -36,7 +71,6 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): cdef double tmp for i in xrange(max_iter): changes = 0 - print("iteration %d" % i) current_mean = means.data # assign pixels to means for k in xrange(n_means): @@ -65,6 +99,7 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): 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()) From 026b6b1df06f0af452b5ed33b3a6d3cfffbd943f Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 3 Aug 2012 11:43:03 +0100 Subject: [PATCH 38/63] Fix initialization in km_segmentation, prettier examples --- doc/examples/plot_km_segmentation.py | 25 ++++-------------- doc/examples/plot_quickshift.py | 31 ++++++++--------------- skimage/__init__.py | 1 + skimage/filter/tests/test_thresholding.py | 2 ++ skimage/segmentation/__init__.py | 3 ++- skimage/segmentation/boundaries.py | 19 ++++++++++++++ skimage/segmentation/km_segmentation.pyx | 11 +++++--- 7 files changed, 46 insertions(+), 46 deletions(-) create mode 100644 skimage/segmentation/boundaries.py diff --git a/doc/examples/plot_km_segmentation.py b/doc/examples/plot_km_segmentation.py index ed8cf983..3a5d6128 100644 --- a/doc/examples/plot_km_segmentation.py +++ b/doc/examples/plot_km_segmentation.py @@ -6,31 +6,16 @@ import matplotlib.pyplot as plt import numpy as np from skimage.data import lena -from skimage.segmentation import km_segmentation +from skimage.segmentation import km_segmentation, visualize_boundaries from skimage.util import img_as_float +from skimage.color import rgb2lab img = img_as_float(lena()).copy("C") -segments = km_segmentation(img, ratio=2.0, n_segments=200) +segments = km_segmentation(rgb2lab(img), ratio=10.0, n_segments=1000) print("number of segments: %d" % len(np.unique(segments))) -plt.subplot(131, title="original") -plt.imshow(img, interpolation='nearest') +boundaries_mine = visualize_boundaries(img, segments) +plt.imshow(boundaries_mine) plt.axis("off") - -plt.subplot(132, title="superpixels") -# shuffle the labels for better visualization -plt.imshow(segments, interpolation='nearest', cmap=plt.cm.prism) -plt.axis("off") - -plt.subplot(133, title="mean color") -colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in - xrange(img.shape[2])] -counts = np.bincount(segments.ravel()) -colors = np.vstack(colors) / counts -plt.imshow(colors.T[segments], interpolation='nearest') -plt.axis("off") - -plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, - bottom=0.02, left=0.02, right=0.98) plt.show() diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py index 5adf1969..9b096459 100644 --- a/doc/examples/plot_quickshift.py +++ b/doc/examples/plot_quickshift.py @@ -25,30 +25,19 @@ import matplotlib.pyplot as plt import numpy as np from skimage.data import lena -from skimage.segmentation import quickshift +from skimage.segmentation import quickshift, visualize_boundaries from skimage.util import img_as_float +from skimage.color import rgb2lab img = img_as_float(lena())[::2, ::2, :].copy("C") -segments = quickshift(img, kernel_size=5, max_dist=20) +segments = quickshift(rgb2lab(img), kernel_size=5, max_dist=20) +segments_rgb = quickshift(img, kernel_size=5, max_dist=20) print("number of segments: %d" % len(np.unique(segments))) - -fig, (ax_org, ax_sp, ax_mean) = plt.subplots(1, 3) -ax_org.set_title("original") -ax_org.imshow(img, interpolation='nearest') -ax_org.axis("off") - -ax_sp.set_title("superpixels") -ax_sp.imshow(segments, interpolation='nearest', cmap=plt.cm.prism) -ax_sp.axis("off") - -colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in - xrange(img.shape[2])] -counts = np.bincount(segments.ravel()) -colors = np.vstack(colors) / counts -ax_mean.set_title("mean color") -ax_mean.imshow(colors.T[segments], interpolation='nearest') -ax_mean.axis("off") -fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, - bottom=0.02, left=0.02, right=0.98) +boundaries = visualize_boundaries(img, segments) +boundaries_rgb = visualize_boundaries(img, segments_rgb) +plt.imshow(boundaries) +plt.figure() +plt.imshow(boundaries_rgb) +plt.axis("off") plt.show() diff --git a/skimage/__init__.py b/skimage/__init__.py index ad37f425..b1c331ff 100644 --- a/skimage/__init__.py +++ b/skimage/__init__.py @@ -88,6 +88,7 @@ test = _setup_test() test_verbose = _setup_test(verbose=True) + def get_log(name=None): """Return a console logger. diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 97d3d9e3..d17f9b84 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -77,11 +77,13 @@ def test_otsu_camera_image(): assert 86 < threshold_otsu(camera) < 88 + def test_otsu_coins_image(): coins = skimage.img_as_ubyte(data.coins()) assert 106 < threshold_otsu(coins) < 108 + def test_otsu_coins_image_as_float(): coins = skimage.img_as_float(data.coins()) assert 0.41 < threshold_otsu(coins) < 0.42 diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index b2c97448..380e937d 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -2,6 +2,7 @@ from .random_walker_segmentation import random_walker from .felzenszwalb import felzenszwalb_segmentation from .km_segmentation import km_segmentation from .quickshift import quickshift +from .boundaries import find_boundaries, visualize_boundaries __all__ = [random_walker, quickshift, felzenszwalb_segmentation, - km_segmentation] + km_segmentation, find_boundaries, visualize_boundaries] diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py new file mode 100644 index 00000000..b40ecb01 --- /dev/null +++ b/skimage/segmentation/boundaries.py @@ -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 diff --git a/skimage/segmentation/km_segmentation.pyx b/skimage/segmentation/km_segmentation.pyx index d0674be4..47faab7f 100644 --- a/skimage/segmentation/km_segmentation.pyx +++ b/skimage/segmentation/km_segmentation.pyx @@ -5,12 +5,12 @@ from scipy import ndimage from ..util import img_as_float -def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): +def km_segmentation(image, n_segments=100, ratio=10., max_iter=10, sigma=1): """Segments image using k-means clustering in Color-(x,y) space. Parameters ---------- - image: (width, height, 3) ndarray + image: (width, height, 3) ndarray Input image ratio: float Balances color-space proximity and image-space proximity. @@ -50,7 +50,8 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): means_y = grid_y[::step, ::step] means_x = grid_x[::step, ::step] - means_color = image[means_y, means_x, :] + n_seeds = len(means_y) + means_color = np.zeros((n_seeds, n_seeds, 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 @@ -63,13 +64,14 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): 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.ones((height, width), dtype=np.float) * np.inf + cdef np.ndarray[dtype=np.float_t, ndim=2] distance = np.empty((height, width)) cdef np.float_t* image_p = image_yx.data cdef np.float_t* distance_p = distance.data cdef np.float_t* current_distance cdef np.float_t* current_pixel cdef double tmp for i in xrange(max_iter): + distance.fill(np.inf) changes = 0 current_mean = means.data # assign pixels to means @@ -105,5 +107,6 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=100, sigma=1.0): means_list = [np.bincount(nearest_mean.ravel(), image_yx[:, :, j].ravel()) for j in xrange(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 From cd1007a0bc3ad88def99ad58acb709e72c341490 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 3 Aug 2012 11:57:02 +0100 Subject: [PATCH 39/63] Rename km_segmentation to slic. They have a PAMI paper now so I guess we should use their name. --- doc/examples/{plot_km_segmentation.py => plot_slic.py} | 4 ++-- skimage/segmentation/__init__.py | 4 ++-- skimage/segmentation/setup.py | 4 ++-- skimage/segmentation/{km_segmentation.pyx => slic.pyx} | 8 ++++---- 4 files changed, 10 insertions(+), 10 deletions(-) rename doc/examples/{plot_km_segmentation.py => plot_slic.py} (73%) rename skimage/segmentation/{km_segmentation.pyx => slic.pyx} (93%) diff --git a/doc/examples/plot_km_segmentation.py b/doc/examples/plot_slic.py similarity index 73% rename from doc/examples/plot_km_segmentation.py rename to doc/examples/plot_slic.py index 3a5d6128..aa743276 100644 --- a/doc/examples/plot_km_segmentation.py +++ b/doc/examples/plot_slic.py @@ -6,12 +6,12 @@ import matplotlib.pyplot as plt import numpy as np from skimage.data import lena -from skimage.segmentation import km_segmentation, visualize_boundaries +from skimage.segmentation import slic, visualize_boundaries from skimage.util import img_as_float from skimage.color import rgb2lab img = img_as_float(lena()).copy("C") -segments = km_segmentation(rgb2lab(img), ratio=10.0, n_segments=1000) +segments = slic(rgb2lab(img), ratio=10.0, n_segments=1000) print("number of segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 380e937d..2fb6902c 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,8 +1,8 @@ from .random_walker_segmentation import random_walker from .felzenszwalb import felzenszwalb_segmentation -from .km_segmentation import km_segmentation +from .slic import slic from .quickshift import quickshift from .boundaries import find_boundaries, visualize_boundaries __all__ = [random_walker, quickshift, felzenszwalb_segmentation, - km_segmentation, find_boundaries, visualize_boundaries] + slic, find_boundaries, visualize_boundaries] diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index 0be6b748..7b7d9cf1 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -17,8 +17,8 @@ def configuration(parent_package='', top_path=None): cython(['quickshift.pyx'], working_path=base_path) config.add_extension('quickshift', sources=['quickshift.c'], include_dirs=[get_numpy_include_dirs()]) - cython(['km_segmentation.pyx'], working_path=base_path) - config.add_extension('km_segmentation', sources=['km_segmentation.c'], + cython(['slic.pyx'], working_path=base_path) + config.add_extension('slic', sources=['slic.c'], include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/segmentation/km_segmentation.pyx b/skimage/segmentation/slic.pyx similarity index 93% rename from skimage/segmentation/km_segmentation.pyx rename to skimage/segmentation/slic.pyx index 47faab7f..3fcd094a 100644 --- a/skimage/segmentation/km_segmentation.pyx +++ b/skimage/segmentation/slic.pyx @@ -5,7 +5,7 @@ from scipy import ndimage from ..util import img_as_float -def km_segmentation(image, n_segments=100, ratio=10., max_iter=10, sigma=1): +def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1): """Segments image using k-means clustering in Color-(x,y) space. Parameters @@ -32,9 +32,9 @@ def km_segmentation(image, n_segments=100, ratio=10., max_iter=10, sigma=1): References ---------- - .. [1] Slic superpixels, Achanta, R. and Shaji, A. and Smith, K. and Lucchi, - A. and Fua, P. and Suesstrunk, S. - Technical Report 2010 + .. [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) From b6059b5672ccff80c9e70d58aa010ef4ac55a076 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 3 Aug 2012 12:00:39 +0100 Subject: [PATCH 40/63] Put RGB2Lab into slic as it seems to be essential. --- doc/examples/plot_slic.py | 3 +-- skimage/segmentation/slic.pyx | 9 +++++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/doc/examples/plot_slic.py b/doc/examples/plot_slic.py index aa743276..8ff59d98 100644 --- a/doc/examples/plot_slic.py +++ b/doc/examples/plot_slic.py @@ -8,10 +8,9 @@ import numpy as np from skimage.data import lena from skimage.segmentation import slic, visualize_boundaries from skimage.util import img_as_float -from skimage.color import rgb2lab img = img_as_float(lena()).copy("C") -segments = slic(rgb2lab(img), ratio=10.0, n_segments=1000) +segments = slic(img, ratio=10.0, n_segments=1000) print("number of segments: %d" % len(np.unique(segments))) diff --git a/skimage/segmentation/slic.pyx b/skimage/segmentation/slic.pyx index 3fcd094a..1430ba63 100644 --- a/skimage/segmentation/slic.pyx +++ b/skimage/segmentation/slic.pyx @@ -3,9 +3,10 @@ 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): +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 @@ -19,6 +20,9 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1): maximum number of iterations of k-means sigma: float Width of Gaussian smoothing kernel for preprocessing. + 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 ------- @@ -28,7 +32,6 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1): Notes ----- The image is smoothed using a Gaussian kernel prior to segmentation. - Best results are achieved if the image is given in Lab color space. References ---------- @@ -41,6 +44,8 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1): if image.shape[2] != 3: ValueError("Only 3-channel 2d images are supported.") image = ndimage.gaussian_filter(img_as_float(image), sigma) + if convert2lab: + image = rgb2lab(image) # initialize on grid: height, width = image.shape[:2] From 1746d2c38c0c059c4d87879a3bcc7d92bc65ed98 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 4 Aug 2012 18:03:03 +0100 Subject: [PATCH 41/63] Some scaling issues to be closer to reference implementation --- skimage/segmentation/quickshift.pyx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index b0a9b2d6..e37463b3 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -65,7 +65,6 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r # we get crazy memory overhead (width * height * windowsize**2) # TODO do smoothing beforehand? - # TODO manage borders somehow? # TODO join orphant roots? # window size for neighboring pixels to consider @@ -92,7 +91,7 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r for c in xrange(channels): dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 dist += (x - x_)**2 + (y - y_)**2 - densities[x, y] += exp(-dist / kernel_size) + densities[x, y] += exp(-dist / (2 * kernel_size**2)) current_pixel_p += channels # this will break ties that otherwise would give us headache @@ -117,7 +116,7 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r if dist < closest: closest = dist parent[x, y] = x_ * width + y_ - dist_parent[x, y] = closest + dist_parent[x, y] = np.sqrt(closest) current_pixel_p += channels dist_parent_flat = dist_parent.ravel() @@ -130,5 +129,5 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r flat = np.unique(flat, return_inverse=True)[1] flat = flat.reshape(width, height) if return_tree: - return flat, parent + return flat, parent, dist_parent return flat From d770fc714d0bb43ac47a7c3e36b02c9791210efe Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 4 Aug 2012 19:56:22 +0100 Subject: [PATCH 42/63] quickshift: convert to lab in function, some comments in code. --- skimage/segmentation/quickshift.pyx | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index e37463b3..3fcc94f7 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -5,6 +5,7 @@ cimport cython from itertools import product from ..util import img_as_float +from ..color import rgb2lab cdef extern from "math.h": @@ -14,7 +15,7 @@ cdef extern from "math.h": @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, random_seed=None): +def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, 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. @@ -34,6 +35,9 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r Higher means less clusters. return_tree: bool Whether to return the full segmentation hierarchy tree + 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 @@ -44,7 +48,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r Notes ----- - The authors advocate to convert the image to Lab color space prior to segmentation. + 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 ---------- @@ -53,8 +58,13 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r """ - image = np.atleast_3d(image) - cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = np.ascontiguousarray(img_as_float(image)) * ratio + 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) + + cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c = np.ascontiguousarray(image) * ratio if random_seed is None: random_state = np.random.RandomState() @@ -64,13 +74,16 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r # We compute the distances twice since otherwise # we get crazy memory overhead (width * height * windowsize**2) - # TODO do smoothing beforehand? # TODO join orphant 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(2 * kernel_size) + cdef int w = int(3 * kernel_size) cdef int width = image_c.shape[0] cdef int height = image_c.shape[1] @@ -95,8 +108,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r current_pixel_p += channels # this will break ties that otherwise would give us headache - densities += random_state.normal(scale=0.00001, size=(width, height)) + # default parent to self: cdef np.ndarray[dtype=np.int_t, ndim=2] parent = np.arange(width * height).reshape(width, height) cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) @@ -121,8 +134,10 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, r dist_parent_flat = dist_parent.ravel() flat = parent.ravel() + # remove parents with distance > max_dist flat[dist_parent_flat > max_dist] = np.arange(width * height)[dist_parent_flat > max_dist] old = np.zeros_like(flat) + # flatten forest (mark each pixel with root of corresponding tree) while (old != flat).any(): old = flat flat = flat[flat] From 10934cfaff34e427481871c2ffa43939cb714896 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 4 Aug 2012 20:27:18 +0100 Subject: [PATCH 43/63] Fixed bug in SLIC smoothing --- skimage/segmentation/slic.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/slic.pyx b/skimage/segmentation/slic.pyx index 1430ba63..37c6b79c 100644 --- a/skimage/segmentation/slic.pyx +++ b/skimage/segmentation/slic.pyx @@ -19,7 +19,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, convert2lab=Tru max_iter: int maximum number of iterations of k-means sigma: float - Width of Gaussian smoothing kernel for preprocessing. + 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. @@ -43,7 +43,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, convert2lab=Tru 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) + image = ndimage.gaussian_filter(img_as_float(image), [sigma, sigma, 0]) if convert2lab: image = rgb2lab(image) From e65d6f6635b5d8289e2b5e5ac6565ea2fa3553bc Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sat, 4 Aug 2012 20:27:31 +0100 Subject: [PATCH 44/63] Added smoothing option to quickshift --- skimage/segmentation/quickshift.pyx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 3fcc94f7..be805a37 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -3,6 +3,7 @@ cimport numpy as np cimport cython from itertools import product +from scipy import ndimage from ..util import img_as_float from ..color import rgb2lab @@ -15,7 +16,7 @@ cdef extern from "math.h": @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, convert2lab=True, random_seed=None): +def quickshift(image, ratio=1., 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. @@ -34,7 +35,9 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, c Cut-off point for data distances. Higher means less clusters. return_tree: bool - Whether to return the full segmentation hierarchy tree + 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. @@ -64,6 +67,7 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, c 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 if random_seed is None: From ec6a1769faa51a8f6977c657e137c6c8d028f804 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 13:51:55 +0100 Subject: [PATCH 45/63] merge segmentation examples --- .../plot_felzenszwalb_segmentation.py | 51 ----------- doc/examples/plot_quickshift.py | 43 ---------- doc/examples/plot_segmentations.py | 85 +++++++++++++++++++ doc/examples/plot_slic.py | 20 ----- 4 files changed, 85 insertions(+), 114 deletions(-) delete mode 100644 doc/examples/plot_felzenszwalb_segmentation.py delete mode 100644 doc/examples/plot_quickshift.py create mode 100644 doc/examples/plot_segmentations.py delete mode 100644 doc/examples/plot_slic.py diff --git a/doc/examples/plot_felzenszwalb_segmentation.py b/doc/examples/plot_felzenszwalb_segmentation.py deleted file mode 100644 index 18ac0f80..00000000 --- a/doc/examples/plot_felzenszwalb_segmentation.py +++ /dev/null @@ -1,51 +0,0 @@ -""" -================================================= -Felzenszwalb's efficient graph based segmentation -================================================= - -This fast 2d image segmentation algorithm, proposed in [1]_ is popular in the -computer vision community. It is often used to extract "superpixels", small -homogeneous image regions, which build the basis for further processing. - -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 -""" -print __doc__ - -import matplotlib.pyplot as plt -import numpy as np - -from skimage.data import lena -from skimage.segmentation import felzenszwalb_segmentation -from skimage.util import img_as_float - -img = img_as_float(lena()) -segments = felzenszwalb_segmentation(img, scale=1) -segments = np.unique(segments, return_inverse=True)[1].reshape(img.shape[:2]) - -print("number of segments: %d" % len(np.unique(segments))) - - -fig, (ax_org, ax_sp, ax_mean) = plt.subplots(1, 3) -ax_org.set_title("original") -ax_org.imshow(img, interpolation='nearest') -ax_org.axis("off") - -ax_sp.set_title("superpixels") -ax_sp.imshow(segments, interpolation='nearest', cmap=plt.cm.prism) -ax_sp.axis("off") - -colors = [np.bincount(segments.ravel(), img[:, :, c].ravel()) for c in - xrange(img.shape[2])] -counts = np.bincount(segments.ravel()) -colors = np.vstack(colors) / counts -ax_mean.set_title("mean color") -ax_mean.imshow(colors.T[segments], interpolation='nearest') -ax_mean.axis("off") -fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, - bottom=0.02, left=0.02, right=0.98) -plt.show() diff --git a/doc/examples/plot_quickshift.py b/doc/examples/plot_quickshift.py deleted file mode 100644 index 9b096459..00000000 --- a/doc/examples/plot_quickshift.py +++ /dev/null @@ -1,43 +0,0 @@ -""" -============================= -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 color+coordinate space, -see [1]_ It is often used to extract "superpixels", small homogeneous image -regions, which build the basis for further processing. - -One of the benefits of quickshift is that it actually computes a -hierarchical segmentation on multiple scales simultaneously. - -Quickshift has two parameters, one controlling the scale of the local -density approximation, the other selecting a level in the hierarchical -segmentation that is produced. - -.. [1] Quick shift and kernel methods for mode seeking, Vedaldi, A. and Soatto, S. - European Conference on Computer Vision, 2008 -""" -print __doc__ - -import matplotlib.pyplot as plt -import numpy as np - -from skimage.data import lena -from skimage.segmentation import quickshift, visualize_boundaries -from skimage.util import img_as_float -from skimage.color import rgb2lab - -img = img_as_float(lena())[::2, ::2, :].copy("C") -segments = quickshift(rgb2lab(img), kernel_size=5, max_dist=20) -segments_rgb = quickshift(img, kernel_size=5, max_dist=20) - -print("number of segments: %d" % len(np.unique(segments))) -boundaries = visualize_boundaries(img, segments) -boundaries_rgb = visualize_boundaries(img, segments_rgb) -plt.imshow(boundaries) -plt.figure() -plt.imshow(boundaries_rgb) -plt.axis("off") -plt.show() diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py new file mode 100644 index 00000000..b204290e --- /dev/null +++ b/doc/examples/plot_segmentations.py @@ -0,0 +1,85 @@ +""" +==================================================== +Comparison of segmentation and superpixel algorithms +==================================================== + +This example compares three popular low-level image segmentation methods. +As it is difficult do obtain good segmentations, and the definition of "good" +often depends on the application, these methods are usually used +for optaining an oversegmentation, also known as superpixels. These superpixels +then serve as the level of operation 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 color+coordinate space, +see [2]_. + +One of the benefits of quickshift is that it actually computes a +hierarchical segmentation on multiple scales simultaneously. + +Quickshift has two parameters, one controlling the scale of the local +density approximation, the other selecting a level in the hierarchical +segmentation that is produced. + +.. [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-kmeans in the 5d color-coordinate space 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. + +.. [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. +""" +print __doc__ + +import matplotlib.pyplot as plt +import numpy as np + +from skimage.data import lena +from skimage.segmentation import felzenszwalb_segmentation, \ + visualize_boundaries, slic, quickshift +from skimage.util import img_as_float + +img = img_as_float(lena()[::2, ::2]) +segments_fz = felzenszwalb_segmentation(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) + +ax[0].imshow(visualize_boundaries(img, segments_fz)) +ax[1].imshow(visualize_boundaries(img, segments_slic)) +ax[2].imshow(visualize_boundaries(img, segments_quick)) +for a in ax: + a.set_xticks(()) + a.set_yticks(()) +plt.show() diff --git a/doc/examples/plot_slic.py b/doc/examples/plot_slic.py deleted file mode 100644 index 8ff59d98..00000000 --- a/doc/examples/plot_slic.py +++ /dev/null @@ -1,20 +0,0 @@ -""" -""" -print __doc__ - -import matplotlib.pyplot as plt -import numpy as np - -from skimage.data import lena -from skimage.segmentation import slic, visualize_boundaries -from skimage.util import img_as_float - -img = img_as_float(lena()).copy("C") -segments = slic(img, ratio=10.0, n_segments=1000) - -print("number of segments: %d" % len(np.unique(segments))) - -boundaries_mine = visualize_boundaries(img, segments) -plt.imshow(boundaries_mine) -plt.axis("off") -plt.show() From b1eb4265ed5704d8c77d83e9972a6c042879df76 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 13:52:15 +0100 Subject: [PATCH 46/63] Post-processing for felzenszwalbs algorithm. --- skimage/segmentation/_felzenszwalb.pyx | 16 ++++++++++++++-- skimage/segmentation/felzenszwalb.py | 7 +++++-- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index a677020a..b66fbd16 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -7,7 +7,7 @@ from skimage.morphology.ccomp cimport find_root, join_trees from ..util import img_as_float -def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): +def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): """Computes Felsenszwalb's efficient graph based segmentation for a single channel. Produces an oversegmentation of a 2d image using a fast, minimum spanning @@ -28,6 +28,8 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): 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 ------- @@ -38,7 +40,8 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): raise ValueError("This algorithm works only on single-channel 2d images." "Got image of shape %s" % str(image.shape)) image = img_as_float(image) - scale = float(scale) + # 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: @@ -88,6 +91,15 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8): segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] cint[seg_new] = costs_p[0] + # postprocessing to remove small segments + edges_p = 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) diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py index cfcf4f23..7da5863a 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.py @@ -4,7 +4,7 @@ import numpy as np from ._felzenszwalb import _felzenszwalb_segmentation_grey -def felzenszwalb_segmentation(image, scale=1, sigma=0.8): +def felzenszwalb_segmentation(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 @@ -29,6 +29,8 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8): 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 ------- @@ -59,7 +61,8 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8): # compute quickshift for each channel for c in xrange(n_channels): channel = np.ascontiguousarray(image[:, :, c]) - s = _felzenszwalb_segmentation_grey(channel, scale=scale, sigma=sigma) + s = _felzenszwalb_segmentation_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 From f22ae2871a4b50dbab76f2ea99cebb4b25d3df6c Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 13:54:02 +0100 Subject: [PATCH 47/63] Pep8 --- doc/examples/plot_segmentations.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index b204290e..9ae45d60 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -3,12 +3,11 @@ Comparison of segmentation and superpixel algorithms ==================================================== -This example compares three popular low-level image segmentation methods. -As it is difficult do obtain good segmentations, and the definition of "good" -often depends on the application, these methods are usually used -for optaining an oversegmentation, also known as superpixels. These superpixels -then serve as the level of operation for more sophisticated algorithms such as CRFs. - +This example compares three popular low-level image segmentation methods. As +it is difficult do obtain good segmentations, and the definition of "good" +often depends on the application, these methods are usually used for optaining +an oversegmentation, also known as superpixels. These superpixels then serve as +the level of operation for more sophisticated algorithms such as CRFs. Felzenszwalb's efficient graph based segmentation From 809871803697e721779be24bad84753356c0838e Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 14:33:30 +0100 Subject: [PATCH 48/63] Fixup tests, add test for slic. --- skimage/segmentation/tests/test_quickshift.py | 19 +++++++------- skimage/segmentation/tests/test_slic.py | 26 +++++++++++++++++++ 2 files changed, 36 insertions(+), 9 deletions(-) create mode 100644 skimage/segmentation/tests/test_slic.py diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py index a904837c..b4bc1e86 100644 --- a/skimage/segmentation/tests/test_quickshift.py +++ b/skimage/segmentation/tests/test_quickshift.py @@ -7,17 +7,17 @@ from skimage.segmentation import quickshift def test_grey(): rnd = np.random.RandomState(0) img = np.zeros((20, 20)) - img[:10, :10] = 0.2 + 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, random_seed=0) + 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], 40) + assert_greater(hist[i], 20) def test_color(): @@ -26,20 +26,21 @@ def test_color(): img[:10, :10, 0] = 1 img[10:, :10, 1] = 1 img[10:, 10:, 2] = 1 - img += 0.2 * rnd.normal(size=img.shape) + img += 0.01 * rnd.normal(size=img.shape) img[img > 1] = 1 img[img < 0] = 0 - seg = quickshift(img, random_seed=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], 2) assert_array_equal(seg[:10, 10:], 1) - assert_array_equal(seg[10:, 10:], 2) + assert_array_equal(seg[10:, 10:], 3) - seg2 = quickshift(img, kernel_size=1, max_dist=3, random_seed=0) + seg2 = quickshift(img, kernel_size=1, max_dist=2, random_seed=0, + convert2lab=False, sigma=0) # very oversegmented: - assert_equal(len(np.unique(seg2)), 18) + assert_equal(len(np.unique(seg2)), 11) # still don't cross lines assert_true((seg2[9, :] != seg2[10, :]).all()) assert_true((seg2[:, 9] != seg2[:, 10]).all()) diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py new file mode 100644 index 00000000..b4d4233a --- /dev/null +++ b/skimage/segmentation/tests/test_slic.py @@ -0,0 +1,26 @@ +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, 20, 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: + 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() From 1db4ebcecb3052386f2d4d30174816a313a94c82 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 20:25:26 +0100 Subject: [PATCH 49/63] MISC some typos in Example, titles set. --- doc/examples/plot_segmentations.py | 34 ++++++++++++++++++------------ 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 9ae45d60..0435d7e5 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -4,10 +4,10 @@ Comparison of segmentation and superpixel algorithms ==================================================== This example compares three popular low-level image segmentation methods. As -it is difficult do obtain good segmentations, and the definition of "good" -often depends on the application, these methods are usually used for optaining +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 -the level of operation for more sophisticated algorithms such as CRFs. +a basis for more sophisticated algorithms such as CRFs. Felzenszwalb's efficient graph based segmentation @@ -26,16 +26,17 @@ 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 color+coordinate space, -see [2]_. +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. see [2]_. One of the benefits of quickshift is that it actually computes a hierarchical segmentation on multiple scales simultaneously. -Quickshift has two parameters, one controlling the scale of the local -density approximation, the other selecting a level in the hierarchical -segmentation that is produced. +Quickshift has three parameters: ``sigma`` controls the scale of the local +density approximation, ``max_dist`` other selecting 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. @@ -44,11 +45,13 @@ segmentation that is produced. SLIC - K-Means based image segmentation --------------------------------------- -This algorithm simply performs K-kmeans in the 5d color-coordinate space 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. +This algorithm simply performs K-kmeans 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 @@ -76,8 +79,11 @@ print("Quickshift number of segments: %d" % len(np.unique(segments_quick))) fig, ax = plt.subplots(1, 3) 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(()) From 530f2c31c402a0609cad80481e5d9ef258ac2599 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 20:25:34 +0100 Subject: [PATCH 50/63] pep8 --- skimage/__init__.py | 1 - skimage/filter/tests/test_thresholding.py | 2 -- skimage/segmentation/__init__.py | 3 -- skimage/segmentation/_felzenszwalb.pyx | 20 +++++++------ skimage/segmentation/quickshift.pyx | 35 ++++++++++++++--------- skimage/segmentation/slic.pyx | 31 ++++++++++++-------- 6 files changed, 53 insertions(+), 39 deletions(-) diff --git a/skimage/__init__.py b/skimage/__init__.py index b1c331ff..ad37f425 100644 --- a/skimage/__init__.py +++ b/skimage/__init__.py @@ -88,7 +88,6 @@ test = _setup_test() test_verbose = _setup_test(verbose=True) - def get_log(name=None): """Return a console logger. diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index d17f9b84..97d3d9e3 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -77,13 +77,11 @@ def test_otsu_camera_image(): assert 86 < threshold_otsu(camera) < 88 - def test_otsu_coins_image(): coins = skimage.img_as_ubyte(data.coins()) assert 106 < threshold_otsu(coins) < 108 - def test_otsu_coins_image_as_float(): coins = skimage.img_as_float(data.coins()) assert 0.41 < threshold_otsu(coins) < 0.42 diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 2fb6902c..b1e6f783 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -3,6 +3,3 @@ from .felzenszwalb import felzenszwalb_segmentation from .slic import slic from .quickshift import quickshift from .boundaries import find_boundaries, visualize_boundaries - -__all__ = [random_walker, quickshift, felzenszwalb_segmentation, - slic, find_boundaries, visualize_boundaries] diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index b66fbd16..f4069274 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -8,7 +8,7 @@ from ..util import img_as_float def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): - """Computes Felsenszwalb's efficient graph based segmentation for a single channel. + """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 parameter ``scale`` sets an @@ -37,8 +37,8 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): 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)) + 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. @@ -49,16 +49,19 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): 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) + 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: width, height = image.shape[:2] - cdef np.ndarray[np.int_t, ndim=2] segments = np.arange(width * height).reshape(width, height) + cdef np.ndarray[np.int_t, ndim=2] segments \ + = np.arange(width * height).reshape(width, height) 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]) + 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) @@ -67,7 +70,8 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): cdef np.int_t *segments_p = segments.data cdef np.int_t *edges_p = edges.data cdef np.float_t *costs_p = costs.data - cdef np.ndarray[np.int_t, ndim=1] segment_size = np.ones(width * height, dtype=np.int) + 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 diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index be805a37..37c9294b 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -16,14 +16,16 @@ cdef extern from "math.h": @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, sigma=0, convert2lab=True, random_seed=None): +def quickshift(image, ratio=1., 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. + Produces an oversegmentation of the image using the quickshift mode-seeking + algorithm. Parameters ---------- - image: (width, height, channels) ndarray + image: (width, height, channels) ndarray Input image ratio: float, between 0 and 1. Balances color-space proximity and image-space proximity. @@ -39,8 +41,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, s 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. + 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 @@ -51,12 +53,14 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, s 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. + 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. + .. [1] Quick shift and kernel methods for mode seeking, + Vedaldi, A. and Soatto, S. European Conference on Computer Vision, 2008 @@ -68,7 +72,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, s 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 + cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c \ + = np.ascontiguousarray(image) * ratio if random_seed is None: random_state = np.random.RandomState() @@ -98,7 +103,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, s cdef np.float_t* image_p = image_c.data cdef np.float_t* current_pixel_p = image_p - cdef np.ndarray[dtype=np.float_t, ndim=2] densities = np.zeros((width, height)) + cdef np.ndarray[dtype=np.float_t, ndim=2] densities \ + = np.zeros((width, height)) # compute densities for x, y in product(xrange(width), xrange(height)): x_min, x_max = max(x - w, 0), min(x + w + 1, width) @@ -115,8 +121,10 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, s densities += random_state.normal(scale=0.00001, size=(width, height)) # default parent to self: - cdef np.ndarray[dtype=np.int_t, ndim=2] parent = np.arange(width * height).reshape(width, height) - cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent = np.zeros((width, height)) + cdef np.ndarray[dtype=np.int_t, ndim=2] parent \ + = np.arange(width * height).reshape(width, height) + cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent \ + = np.zeros((width, height)) # find nearest node with higher density current_pixel_p = image_p for x, y in product(xrange(width), xrange(height)): @@ -139,7 +147,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, s dist_parent_flat = dist_parent.ravel() flat = parent.ravel() # remove parents with distance > max_dist - flat[dist_parent_flat > max_dist] = np.arange(width * height)[dist_parent_flat > 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(): diff --git a/skimage/segmentation/slic.pyx b/skimage/segmentation/slic.pyx index 37c6b79c..0d0adc49 100644 --- a/skimage/segmentation/slic.pyx +++ b/skimage/segmentation/slic.pyx @@ -6,7 +6,8 @@ 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): +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 @@ -19,10 +20,12 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, convert2lab=Tru max_iter: int maximum number of iterations of k-means sigma: float - Width of Gaussian smoothing kernel for preprocessing. Zero means no smoothing. + 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. + 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 ------- @@ -57,19 +60,23 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, convert2lab=Tru n_seeds = len(means_y) means_color = np.zeros((n_seeds, n_seeds, 3)) - cdef np.ndarray[dtype=np.float_t, ndim=2] means = np.dstack([means_y, means_x, means_color]).reshape(-1, 5) + 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 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.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 = image_yx.data cdef np.float_t* distance_p = distance.data cdef np.float_t* current_distance @@ -93,8 +100,8 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, convert2lab=Tru mean_entry = current_mean dist_mean = 0 for c in range(5): - # you would think the compiler can optimize this itself. - # mine can't (with O2) + # you would think the compiler can optimize this + # itself. mine can't (with O2) tmp = current_pixel[0] - mean_entry[0] dist_mean += tmp * tmp current_pixel += 1 @@ -109,8 +116,8 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, convert2lab=Tru if changes == 0: break # recompute means: - means_list = [np.bincount(nearest_mean.ravel(), image_yx[:, :, j].ravel()) - for j in xrange(5)] + means_list = [np.bincount(nearest_mean.ravel(), + image_yx[:, :, j].ravel()) for j in xrange(5)] in_mean = np.bincount(nearest_mean.ravel()) in_mean[in_mean == 0] = 1 means = (np.vstack(means_list) / in_mean).T.copy("C") From e0034e384a050cdb8ea2a55fbd10fa60f9b66665 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 20:30:43 +0100 Subject: [PATCH 51/63] Minor fixes addressing @tonysyu's comments. --- skimage/segmentation/_felzenszwalb.pyx | 12 ++++-------- skimage/segmentation/quickshift.pyx | 7 ++----- 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index f4069274..15b528ad 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -11,21 +11,17 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, 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 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. - + 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: (width, height) ndarray + image: ndarray Input image scale: float - Free parameter. Higher means larger clusters. + Sets the obervation level. Higher means larger clusters. sigma: float Width of Gaussian kernel used in preprocessing. min_size: int @@ -33,7 +29,7 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): Returns ------- - segment_mask: ndarray, [width, height] + segment_mask: (height, width) ndarray Integer mask indicating segment labels. """ if image.ndim != 2: diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index 37c9294b..a3c5ec9f 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -75,15 +75,12 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, cdef np.ndarray[dtype=np.float_t, ndim=3, mode="c"] image_c \ = np.ascontiguousarray(image) * ratio - if random_seed is None: - random_state = np.random.RandomState() - else: - random_state = np.random.RandomState(random_seed) + random_state = np.random.RandomState(random_seed) # We compute the distances twice since otherwise # we get crazy memory overhead (width * height * windowsize**2) - # TODO join orphant roots? + # 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 From 73dd46019b4eb7934d4a96a2d23430f07a717b82 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 21:06:37 +0100 Subject: [PATCH 52/63] FIX width/height trouble, add non-regression test --- skimage/segmentation/_felzenszwalb.pyx | 6 ++--- skimage/segmentation/quickshift.pyx | 26 +++++++++---------- skimage/segmentation/slic.pyx | 8 +++--- .../segmentation/tests/test_felzenszwalb.py | 4 +-- skimage/segmentation/tests/test_quickshift.py | 13 +++++----- skimage/segmentation/tests/test_slic.py | 3 ++- 6 files changed, 31 insertions(+), 29 deletions(-) diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index 15b528ad..5c1e097e 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -49,9 +49,9 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): down_cost.ravel(), dright_cost.ravel(), uright_cost.ravel()]).astype(np.float) # compute edges between pixels: - width, height = image.shape[:2] + height, width = image.shape[:2] cdef np.ndarray[np.int_t, ndim=2] segments \ - = np.arange(width * height).reshape(width, height) + = 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()] @@ -107,4 +107,4 @@ def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): old = flat flat = flat[flat] flat = np.unique(flat, return_inverse=True)[1] - return flat.reshape((width, height)) + return flat.reshape((height, width)) diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/quickshift.pyx index a3c5ec9f..dc4d86e1 100644 --- a/skimage/segmentation/quickshift.pyx +++ b/skimage/segmentation/quickshift.pyx @@ -91,8 +91,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, raise ValueError("Sigma should be >= 1") cdef int w = int(3 * kernel_size) - cdef int width = image_c.shape[0] - cdef int height = image_c.shape[1] + cdef int height = image_c.shape[0] + cdef int width = image_c.shape[1] cdef int channels = image_c.shape[2] cdef float closest, dist cdef int x, y, x_, y_ @@ -101,11 +101,11 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, cdef np.float_t* current_pixel_p = image_p cdef np.ndarray[dtype=np.float_t, ndim=2] densities \ - = np.zeros((width, height)) + = np.zeros((height, width)) # compute densities - for x, y in product(xrange(width), xrange(height)): - x_min, x_max = max(x - w, 0), min(x + w + 1, width) - y_min, y_max = max(y - w, 0), min(y + w + 1, height) + for x, y in product(xrange(height), xrange(width)): + x_min, x_max = max(x - w, 0), min(x + w + 1, height) + y_min, y_max = max(y - w, 0), min(y + w + 1, width) for x_, y_ in product(xrange(x_min, x_max), xrange(y_min, y_max)): dist = 0 for c in xrange(channels): @@ -115,20 +115,20 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, current_pixel_p += channels # this will break ties that otherwise would give us headache - densities += random_state.normal(scale=0.00001, size=(width, height)) + 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(width, height) + = np.arange(width * height).reshape(height, width) cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent \ - = np.zeros((width, height)) + = np.zeros((height, width)) # find nearest node with higher density current_pixel_p = image_p - for x, y in product(xrange(width), xrange(height)): + for x, y in product(xrange(height), xrange(width)): current_density = densities[x, y] closest = np.inf - x_min, x_max = max(x - w, 0), min(x + w + 1, width) - y_min, y_max = max(y - w, 0), min(y + w + 1, height) + x_min, x_max = max(x - w, 0), min(x + w + 1, height) + y_min, y_max = max(y - w, 0), min(y + w + 1, width) for x_, y_ in product(xrange(x_min, x_max), xrange(y_min, y_max)): if densities[x_, y_] > current_density: dist = 0 @@ -152,7 +152,7 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, old = flat flat = flat[flat] flat = np.unique(flat, return_inverse=True)[1] - flat = flat.reshape(width, height) + flat = flat.reshape(height, width) if return_tree: return flat, parent, dist_parent return flat diff --git a/skimage/segmentation/slic.pyx b/skimage/segmentation/slic.pyx index 0d0adc49..652f977e 100644 --- a/skimage/segmentation/slic.pyx +++ b/skimage/segmentation/slic.pyx @@ -53,13 +53,13 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, # initialize on grid: height, width = image.shape[:2] # approximate grid size for desired n_segments - step = np.sqrt(height * width / n_segments) + 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] + print(means_y, means_x) - n_seeds = len(means_y) - means_color = np.zeros((n_seeds, n_seeds, 3)) + 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 @@ -92,7 +92,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, 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, height)) + x_max = int(min(current_mean[1] + 2 * step, width)) for y in xrange(y_min, y_max): current_pixel = &image_p[5 * (y * width + x_min)] current_distance = &distance_p[y * width + x_min] diff --git a/skimage/segmentation/tests/test_felzenszwalb.py b/skimage/segmentation/tests/test_felzenszwalb.py index f6cca31b..fe68c443 100644 --- a/skimage/segmentation/tests/test_felzenszwalb.py +++ b/skimage/segmentation/tests/test_felzenszwalb.py @@ -6,7 +6,7 @@ from skimage.segmentation import felzenszwalb_segmentation def test_grey(): # very weak tests. This algorithm is pretty unstable. - img = np.zeros((20, 20)) + img = np.zeros((20, 21)) img[:10, 10:] = 0.2 img[10:, :10] = 0.4 img[10:, 10:] = 0.6 @@ -21,7 +21,7 @@ def test_grey(): def test_color(): # very weak tests. This algorithm is pretty unstable. - img = np.zeros((20, 20, 3)) + img = np.zeros((20, 21, 3)) img[:10, :10, 0] = 1 img[10:, :10, 1] = 1 img[10:, 10:, 2] = 1 diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py index b4bc1e86..5c6eb024 100644 --- a/skimage/segmentation/tests/test_quickshift.py +++ b/skimage/segmentation/tests/test_quickshift.py @@ -6,12 +6,13 @@ from skimage.segmentation import quickshift def test_grey(): rnd = np.random.RandomState(0) - img = np.zeros((20, 20)) + 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) + 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: @@ -22,7 +23,7 @@ def test_grey(): def test_color(): rnd = np.random.RandomState(0) - img = np.zeros((20, 20, 3)) + img = np.zeros((20, 21, 3)) img[:10, :10, 0] = 1 img[10:, :10, 1] = 1 img[10:, 10:, 2] = 1 @@ -33,14 +34,14 @@ def test_color(): # 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], 3) assert_array_equal(seg[:10, 10:], 1) - assert_array_equal(seg[10:, 10:], 3) + 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)), 11) + 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()) diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py index b4d4233a..f2d6698d 100644 --- a/skimage/segmentation/tests/test_slic.py +++ b/skimage/segmentation/tests/test_slic.py @@ -5,7 +5,7 @@ from skimage.segmentation import slic def test_color(): rnd = np.random.RandomState(0) - img = np.zeros((20, 20, 3)) + img = np.zeros((20, 21, 3)) img[:10, :10, 0] = 1 img[10:, :10, 1] = 1 img[10:, 10:, 2] = 1 @@ -14,6 +14,7 @@ def test_color(): 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) From f421587aa4a1e29b891380103f9ddce93aec32d5 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 21:10:29 +0100 Subject: [PATCH 53/63] ENH renamed "felzenszwalb_segmentation" to "felzenszwalb", remove debug output from slic --- doc/examples/plot_segmentations.py | 4 ++-- skimage/segmentation/__init__.py | 2 +- skimage/segmentation/_felzenszwalb.pyx | 2 +- skimage/segmentation/felzenszwalb.py | 8 ++++---- skimage/segmentation/slic.pyx | 1 - 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 0435d7e5..890ffcce 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -63,12 +63,12 @@ import matplotlib.pyplot as plt import numpy as np from skimage.data import lena -from skimage.segmentation import felzenszwalb_segmentation, \ +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_segmentation(img, scale=100, sigma=0.5, min_size=50) +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) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index b1e6f783..cb06108d 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,5 +1,5 @@ from .random_walker_segmentation import random_walker -from .felzenszwalb import felzenszwalb_segmentation +from .felzenszwalb import felzenszwalb from .slic import slic from .quickshift import quickshift from .boundaries import find_boundaries, visualize_boundaries diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/_felzenszwalb.pyx index 5c1e097e..76320d3a 100644 --- a/skimage/segmentation/_felzenszwalb.pyx +++ b/skimage/segmentation/_felzenszwalb.pyx @@ -7,7 +7,7 @@ from skimage.morphology.ccomp cimport find_root, join_trees from ..util import img_as_float -def _felzenszwalb_segmentation_grey(image, scale=1, sigma=0.8, min_size=20): +def _felzenszwalb_grey(image, scale=1, sigma=0.8, 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 diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/felzenszwalb.py index 7da5863a..cf79706d 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/felzenszwalb.py @@ -1,10 +1,10 @@ import warnings import numpy as np -from ._felzenszwalb import _felzenszwalb_segmentation_grey +from ._felzenszwalb import _felzenszwalb_grey -def felzenszwalb_segmentation(image, scale=1, sigma=0.8, min_size=20): +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 @@ -46,7 +46,7 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8, min_size=20): #image = img_as_float(image) if image.ndim == 2: # assume single channel image - return _felzenszwalb_segmentation_grey(image, scale=scale, sigma=sigma) + return _felzenszwalb_grey(image, scale=scale, sigma=sigma) elif image.ndim != 3: raise ValueError("Got image with ndim=%d, don't know" @@ -61,7 +61,7 @@ def felzenszwalb_segmentation(image, scale=1, sigma=0.8, min_size=20): # compute quickshift for each channel for c in xrange(n_channels): channel = np.ascontiguousarray(image[:, :, c]) - s = _felzenszwalb_segmentation_grey(channel, scale=scale, sigma=sigma, + s = _felzenszwalb_grey(channel, scale=scale, sigma=sigma, min_size=min_size) segmentations.append(s) diff --git a/skimage/segmentation/slic.pyx b/skimage/segmentation/slic.pyx index 652f977e..98a6f6bc 100644 --- a/skimage/segmentation/slic.pyx +++ b/skimage/segmentation/slic.pyx @@ -57,7 +57,6 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, grid_y, grid_x = np.mgrid[:height, :width] means_y = grid_y[::step, ::step] means_x = grid_x[::step, ::step] - print(means_y, means_x) means_color = np.zeros((means_y.shape[0], means_y.shape[1], 3)) cdef np.ndarray[dtype=np.float_t, ndim=2] means \ From f56034630999f08a50a15e2687571b39d97294c6 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 21:31:24 +0100 Subject: [PATCH 54/63] Trying to avoid name collisions. --- skimage/segmentation/__init__.py | 6 +++--- .../{felzenszwalb.py => _felzenszwalb.py} | 2 +- .../segmentation/{quickshift.pyx => _quickshift.pyx} | 0 skimage/segmentation/{slic.pyx => _slic.pyx} | 0 .../{_felzenszwalb.pyx => felzenszwalb_cy.pyx} | 0 skimage/segmentation/setup.py | 12 ++++++------ 6 files changed, 10 insertions(+), 10 deletions(-) rename skimage/segmentation/{felzenszwalb.py => _felzenszwalb.py} (98%) rename skimage/segmentation/{quickshift.pyx => _quickshift.pyx} (100%) rename skimage/segmentation/{slic.pyx => _slic.pyx} (100%) rename skimage/segmentation/{_felzenszwalb.pyx => felzenszwalb_cy.pyx} (100%) diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index cb06108d..69a580c6 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,5 +1,5 @@ from .random_walker_segmentation import random_walker -from .felzenszwalb import felzenszwalb -from .slic import slic -from .quickshift import quickshift +from ._felzenszwalb import felzenszwalb +from ._slic import slic +from ._quickshift import quickshift from .boundaries import find_boundaries, visualize_boundaries diff --git a/skimage/segmentation/felzenszwalb.py b/skimage/segmentation/_felzenszwalb.py similarity index 98% rename from skimage/segmentation/felzenszwalb.py rename to skimage/segmentation/_felzenszwalb.py index cf79706d..f84f3e56 100644 --- a/skimage/segmentation/felzenszwalb.py +++ b/skimage/segmentation/_felzenszwalb.py @@ -1,7 +1,7 @@ import warnings import numpy as np -from ._felzenszwalb import _felzenszwalb_grey +from .felzenszwalb_cy import _felzenszwalb_grey def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): diff --git a/skimage/segmentation/quickshift.pyx b/skimage/segmentation/_quickshift.pyx similarity index 100% rename from skimage/segmentation/quickshift.pyx rename to skimage/segmentation/_quickshift.pyx diff --git a/skimage/segmentation/slic.pyx b/skimage/segmentation/_slic.pyx similarity index 100% rename from skimage/segmentation/slic.pyx rename to skimage/segmentation/_slic.pyx diff --git a/skimage/segmentation/_felzenszwalb.pyx b/skimage/segmentation/felzenszwalb_cy.pyx similarity index 100% rename from skimage/segmentation/_felzenszwalb.pyx rename to skimage/segmentation/felzenszwalb_cy.pyx diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index 7b7d9cf1..ec092ffe 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -11,14 +11,14 @@ def configuration(parent_package='', top_path=None): config = Configuration('segmentation', parent_package, top_path) - cython(['_felzenszwalb.pyx'], working_path=base_path) - config.add_extension('_felzenszwalb', sources=['_felzenszwalb.c'], + 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.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'], + cython(['_slic.pyx'], working_path=base_path) + config.add_extension('_slic', sources=['_slic.c'], include_dirs=[get_numpy_include_dirs()]) return config From 75e3067acd3a7f38635b4f8d57855a0b5a6a4849 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 21:39:03 +0100 Subject: [PATCH 55/63] DOC a bit nicer figure --- doc/examples/plot_segmentations.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 890ffcce..6927b323 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -77,6 +77,8 @@ 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(15, 6, 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") From ea02bc6170d96e2bd5b610ae24c73dd77b400efa Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Sun, 5 Aug 2012 23:11:50 +0100 Subject: [PATCH 56/63] make figure smaller again. --- doc/examples/plot_segmentations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 6927b323..8830cebb 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -77,7 +77,7 @@ 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(15, 6, forward=True) +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)) From f88a29b091c6af19e7d3df3bc8f91c49216106e5 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 10 Aug 2012 10:03:20 +0100 Subject: [PATCH 57/63] ENH minor speedups. --- skimage/segmentation/_slic.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index 98a6f6bc..684740d6 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -51,9 +51,10 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, image = rgb2lab(image) # initialize on grid: + cdef int height, width height, width = image.shape[:2] # approximate grid size for desired n_segments - step = np.ceil(np.sqrt(height * width / 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] From 8f5337a2bf82267b14b01630f3e513042b1d8082 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 10 Aug 2012 10:31:23 +0100 Subject: [PATCH 58/63] ENH added some cdefs --- skimage/segmentation/_quickshift.pyx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index dc4d86e1..c050688e 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -16,7 +16,7 @@ cdef extern from "math.h": @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, +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. @@ -94,7 +94,8 @@ def quickshift(image, ratio=1., kernel_size=5, max_dist=10, return_tree=False, cdef int height = image_c.shape[0] cdef int width = image_c.shape[1] cdef int channels = image_c.shape[2] - cdef float closest, dist + cdef double current_density, closest, dist + cdef int x, y, x_, y_ cdef np.float_t* image_p = image_c.data From 312b03d1b141bf6e805544322bcae3f9e34bac86 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 10 Aug 2012 10:35:23 +0100 Subject: [PATCH 59/63] ENH more speeeeed --- skimage/segmentation/_quickshift.pyx | 57 +++++++++++++++------------- 1 file changed, 31 insertions(+), 26 deletions(-) diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index c050688e..4267c17b 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -11,6 +11,7 @@ from ..color import rgb2lab cdef extern from "math.h": double exp(double) + double sqrt(double) @cython.boundscheck(False) @@ -104,16 +105,18 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa cdef np.ndarray[dtype=np.float_t, ndim=2] densities \ = np.zeros((height, width)) # compute densities - for x, y in product(xrange(height), xrange(width)): - x_min, x_max = max(x - w, 0), min(x + w + 1, height) - y_min, y_max = max(y - w, 0), min(y + w + 1, width) - for x_, y_ in product(xrange(x_min, x_max), xrange(y_min, y_max)): - dist = 0 - for c in xrange(channels): - dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 - dist += (x - x_)**2 + (y - y_)**2 - densities[x, y] += exp(-dist / (2 * kernel_size**2)) - current_pixel_p += channels + for x in range(height): + for y in range(width): + x_min, x_max = max(x - w, 0), min(x + w + 1, height) + y_min, y_max = max(y - w, 0), min(y + w + 1, width) + for x_ in range(x_min, x_max): + for y_ in range(y_min, y_max): + dist = 0 + for c in range(channels): + dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 + dist += (x - x_)**2 + (y - y_)**2 + densities[x, y] += 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)) @@ -125,22 +128,24 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa = np.zeros((height, width)) # find nearest node with higher density current_pixel_p = image_p - for x, y in product(xrange(height), xrange(width)): - current_density = densities[x, y] - closest = np.inf - x_min, x_max = max(x - w, 0), min(x + w + 1, height) - y_min, y_max = max(y - w, 0), min(y + w + 1, width) - for x_, y_ in product(xrange(x_min, x_max), xrange(y_min, y_max)): - if densities[x_, y_] > current_density: - dist = 0 - for c in xrange(channels): - dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 - dist += (x - x_)**2 + (y - y_)**2 - if dist < closest: - closest = dist - parent[x, y] = x_ * width + y_ - dist_parent[x, y] = np.sqrt(closest) - current_pixel_p += channels + for x in range(height): + for y in range(width): + current_density = densities[x, y] + closest = np.inf + x_min, x_max = max(x - w, 0), min(x + w + 1, height) + y_min, y_max = max(y - w, 0), min(y + w + 1, width) + for x_ in range(x_min, x_max): + for y_ in range(y_min, y_max): + if densities[x_, y_] > current_density: + dist = 0 + for c in range(channels): + dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 + dist += (x - x_)**2 + (y - y_)**2 + if dist < closest: + closest = dist + parent[x, y] = x_ * width + y_ + dist_parent[x, y] = sqrt(closest) + current_pixel_p += channels dist_parent_flat = dist_parent.ravel() flat = parent.ravel() From 8d769a4cd966a29ad0764add3a73b4c074c003f6 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Fri, 10 Aug 2012 10:48:38 +0100 Subject: [PATCH 60/63] ENH Felzenszwalbs segmentation somewhat faster --- skimage/segmentation/felzenszwalb_cy.pyx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/skimage/segmentation/felzenszwalb_cy.pyx b/skimage/segmentation/felzenszwalb_cy.pyx index 76320d3a..c5f3e705 100644 --- a/skimage/segmentation/felzenszwalb_cy.pyx +++ b/skimage/segmentation/felzenszwalb_cy.pyx @@ -1,13 +1,16 @@ 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 - -def _felzenszwalb_grey(image, scale=1, sigma=0.8, min_size=20): +@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 @@ -70,7 +73,7 @@ def _felzenszwalb_grey(image, scale=1, sigma=0.8, min_size=20): = 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 + 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. From 37c0ffe0724531e6ab27b28038a6e7ba2030f352 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 20 Aug 2012 20:01:41 +0100 Subject: [PATCH 61/63] cosmit: changed the names of x and y to r and c. that was no fun, i tell you. --- skimage/segmentation/_quickshift.pyx | 55 ++++++++++++++-------------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index 4267c17b..a924d72f 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -18,7 +18,7 @@ cdef extern from "math.h": @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): + 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 @@ -78,9 +78,6 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa random_state = np.random.RandomState(random_seed) - # We compute the distances twice since otherwise - # we get crazy memory overhead (width * height * windowsize**2) - # 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. @@ -97,7 +94,7 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa cdef int channels = image_c.shape[2] cdef double current_density, closest, dist - cdef int x, y, x_, y_ + cdef int r, c, r_, c_, channel cdef np.float_t* image_p = image_c.data cdef np.float_t* current_pixel_p = image_p @@ -105,17 +102,17 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa cdef np.ndarray[dtype=np.float_t, ndim=2] densities \ = np.zeros((height, width)) # compute densities - for x in range(height): - for y in range(width): - x_min, x_max = max(x - w, 0), min(x + w + 1, height) - y_min, y_max = max(y - w, 0), min(y + w + 1, width) - for x_ in range(x_min, x_max): - for y_ in range(y_min, y_max): + 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 c in range(channels): - dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 - dist += (x - x_)**2 + (y - y_)**2 - densities[x, y] += exp(-dist / (2 * kernel_size**2)) + 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 @@ -128,23 +125,25 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa = np.zeros((height, width)) # find nearest node with higher density current_pixel_p = image_p - for x in range(height): - for y in range(width): - current_density = densities[x, y] + for r in range(height): + for c in range(width): + current_density = densities[r, c] closest = np.inf - x_min, x_max = max(x - w, 0), min(x + w + 1, height) - y_min, y_max = max(y - w, 0), min(y + w + 1, width) - for x_ in range(x_min, x_max): - for y_ in range(y_min, y_max): - if densities[x_, y_] > current_density: + 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 - for c in range(channels): - dist += (current_pixel_p[c] - image_c[x_, y_, c])**2 - dist += (x - x_)**2 + (y - y_)**2 + # 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[x, y] = x_ * width + y_ - dist_parent[x, y] = sqrt(closest) + parent[r, c] = r_ * width + c_ + dist_parent[r, c] = sqrt(closest) current_pixel_p += channels dist_parent_flat = dist_parent.ravel() From fe2a4334faa0383fde6c29c791508151d632f704 Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 20 Aug 2012 20:22:06 +0100 Subject: [PATCH 62/63] ENH addressed (hopefully all) of Tony's and Stefan's comments. --- doc/examples/plot_segmentations.py | 9 +++--- skimage/segmentation/_felzenszwalb.py | 29 ++++++++++--------- skimage/segmentation/_quickshift.pyx | 28 +++++++++--------- skimage/segmentation/_slic.pyx | 16 +++++----- skimage/segmentation/felzenszwalb_cy.pyx | 7 +++-- .../segmentation/tests/test_felzenszwalb.py | 6 ++-- 6 files changed, 48 insertions(+), 47 deletions(-) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 8830cebb..09012bbd 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -28,13 +28,13 @@ 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. see [2]_. +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 three parameters: ``sigma`` controls the scale of the local -density approximation, ``max_dist`` other selecting a level in the hierarchical +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``. @@ -45,7 +45,7 @@ color-space and distance in image-space, given by ``ratio``. SLIC - K-Means based image segmentation --------------------------------------- -This algorithm simply performs K-kmeans in the 5d space of color information +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 @@ -57,7 +57,6 @@ of Quickshift, while ``n_segments`` chooses the number of centers for kmeans. Pascal Fua, and Sabine Suesstrunk, SLIC Superpixels Compared to State-of-the-art Superpixel Methods, TPAMI, May 2012. """ -print __doc__ import matplotlib.pyplot as plt import numpy as np diff --git a/skimage/segmentation/_felzenszwalb.py b/skimage/segmentation/_felzenszwalb.py index f84f3e56..5729bd95 100644 --- a/skimage/segmentation/_felzenszwalb.py +++ b/skimage/segmentation/_felzenszwalb.py @@ -17,24 +17,24 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): controlled indirectly through ``scale``. Segment size within an image can vary greatly depending on local contrast. - Calls the algorithm on each channel separately, then combines - using "and", i.e. two pixels are in the same segment if they are - in the same segment for each channel. + 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) ndarray - Input image - scale: float + image : (width, height, 3) or (width, height) ndarray + Input image. + scale : float Free parameter. Higher means larger clusters. - sigma: float + sigma : float Width of Gaussian kernel used in preprocessing. - min_size: int + min_size : int Minimum component size. Enforced using postprocessing. Returns ------- - segment_mask: ndarray, [width, height] + segment_mask : (width, height) ndarray Integer mask indicating segment labels. References @@ -49,20 +49,21 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): return _felzenszwalb_grey(image, scale=scale, sigma=sigma) elif image.ndim != 3: - raise ValueError("Got image with ndim=%d, don't know" - " what to do." % image.ndim) + 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]) + " wanted?" % image.shape[2]) segmentations = [] # compute quickshift for each channel for c in xrange(n_channels): channel = np.ascontiguousarray(image[:, :, c]) s = _felzenszwalb_grey(channel, scale=scale, sigma=sigma, - min_size=min_size) + min_size=min_size) segmentations.append(s) # put pixels in same segment only if in the same segment in all images @@ -70,7 +71,7 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): n0 = segmentations[0].max() + 1 n1 = segmentations[1].max() + 1 segmentation = (segmentations[0] + segmentations[1] * n0 - + segmentations[2] * n0 * n1) + + 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]) diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index a924d72f..57009ad1 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -26,30 +26,30 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, return_tree=Fa Parameters ---------- - image: (width, height, channels) ndarray - Input image - ratio: float, between 0 and 1. + 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 + kernel_size : float Width of Gaussian kernel used in smoothing the - sample density. Higher means less clusters. - max_dist: float + sample density. Higher means fewer clusters. + max_dist : float Cut-off point for data distances. - Higher means less clusters. - return_tree: bool + Higher means fewer clusters. + return_tree : bool Whether to return the full segmentation hierarchy tree and distances. - sigma: float + sigma : float Width for Gaussian smoothing as preprocessing. Zero means no smoothing. - convert2lab: bool + 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 + 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: ndarray, [width, height] + segment_mask : (width, height) ndarray Integer mask indicating segment labels. Notes diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index 684740d6..a4f37fb2 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -12,24 +12,24 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, Parameters ---------- - image: (width, height, 3) ndarray - Input image + 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 + 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 + 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: ndarray, [width, height] + segment_mask : (width, height) ndarray Integer mask indicating segment labels. Notes @@ -100,7 +100,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, mean_entry = current_mean dist_mean = 0 for c in range(5): - # you would think the compiler can optimize this + # 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 diff --git a/skimage/segmentation/felzenszwalb_cy.pyx b/skimage/segmentation/felzenszwalb_cy.pyx index c5f3e705..d2c2e00a 100644 --- a/skimage/segmentation/felzenszwalb_cy.pyx +++ b/skimage/segmentation/felzenszwalb_cy.pyx @@ -14,7 +14,7 @@ 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. + 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. @@ -22,11 +22,12 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, int min_size=20): Parameters ---------- image: ndarray - Input image + Input image. scale: float Sets the obervation level. Higher means larger clusters. sigma: float - Width of Gaussian kernel used in preprocessing. + Width of Gaussian smoothing kernel used in preprocessing. + Larger sigma gives smother segment boundaries. min_size: int Minimum component size. Enforced using postprocessing. diff --git a/skimage/segmentation/tests/test_felzenszwalb.py b/skimage/segmentation/tests/test_felzenszwalb.py index fe68c443..ebda2a38 100644 --- a/skimage/segmentation/tests/test_felzenszwalb.py +++ b/skimage/segmentation/tests/test_felzenszwalb.py @@ -1,7 +1,7 @@ import numpy as np from numpy.testing import assert_equal, assert_array_equal from nose.tools import assert_greater -from skimage.segmentation import felzenszwalb_segmentation +from skimage.segmentation import felzenszwalb def test_grey(): @@ -10,7 +10,7 @@ def test_grey(): img[:10, 10:] = 0.2 img[10:, :10] = 0.4 img[10:, 10:] = 0.6 - seg = felzenszwalb_segmentation(img, sigma=0) + seg = felzenszwalb(img, sigma=0) # we expect 4 segments: assert_equal(len(np.unique(seg)), 4) # that mostly respect the 4 regions: @@ -25,7 +25,7 @@ def test_color(): img[:10, :10, 0] = 1 img[10:, :10, 1] = 1 img[10:, 10:, 2] = 1 - seg = felzenszwalb_segmentation(img, sigma=0) + seg = felzenszwalb(img, sigma=0) # we expect 4 segments: assert_equal(len(np.unique(seg)), 4) assert_array_equal(seg[:10, :10], 0) From 6b1dab9f9a7777642b7126b65f2d2ab949d86b5d Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 20 Aug 2012 22:30:58 +0100 Subject: [PATCH 63/63] MISC move felzenszwalb_cy.pyx to _felzenszwalb_cy.pyx, don't use xrange when not necessary --- doc/examples/plot_segmentations.py | 6 +++--- skimage/segmentation/_felzenszwalb.py | 4 ++-- .../{felzenszwalb_cy.pyx => _felzenszwalb_cy.pyx} | 0 skimage/segmentation/_slic.pyx | 12 ++++++------ skimage/segmentation/setup.py | 4 ++-- 5 files changed, 13 insertions(+), 13 deletions(-) rename skimage/segmentation/{felzenszwalb_cy.pyx => _felzenszwalb_cy.pyx} (100%) diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 09012bbd..a8ee7cea 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -12,7 +12,7 @@ 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 +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 @@ -25,9 +25,9 @@ local contrast. Quickshift image segmentation ----------------------------- -Quickshift is a relatively recent 2d image segmentation algorithm, based on an +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 +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 diff --git a/skimage/segmentation/_felzenszwalb.py b/skimage/segmentation/_felzenszwalb.py index 5729bd95..67971a96 100644 --- a/skimage/segmentation/_felzenszwalb.py +++ b/skimage/segmentation/_felzenszwalb.py @@ -1,7 +1,7 @@ import warnings import numpy as np -from .felzenszwalb_cy import _felzenszwalb_grey +from ._felzenszwalb_cy import _felzenszwalb_grey def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): @@ -60,7 +60,7 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): " wanted?" % image.shape[2]) segmentations = [] # compute quickshift for each channel - for c in xrange(n_channels): + for c in range(n_channels): channel = np.ascontiguousarray(image[:, :, c]) s = _felzenszwalb_grey(channel, scale=scale, sigma=sigma, min_size=min_size) diff --git a/skimage/segmentation/felzenszwalb_cy.pyx b/skimage/segmentation/_felzenszwalb_cy.pyx similarity index 100% rename from skimage/segmentation/felzenszwalb_cy.pyx rename to skimage/segmentation/_felzenszwalb_cy.pyx diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index a4f37fb2..ecb58efe 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -45,7 +45,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, """ image = np.atleast_3d(image) if image.shape[2] != 3: - ValueError("Only 3-channel 2d images are supported.") + ValueError("Only 3-channel 2D images are supported.") image = ndimage.gaussian_filter(img_as_float(image), [sigma, sigma, 0]) if convert2lab: image = rgb2lab(image) @@ -82,21 +82,21 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, cdef np.float_t* current_distance cdef np.float_t* current_pixel cdef double tmp - for i in xrange(max_iter): + for i in range(max_iter): distance.fill(np.inf) changes = 0 current_mean = means.data # assign pixels to means - for k in xrange(n_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 xrange(y_min, y_max): + 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 xrange(x_min, x_max): + for x in range(x_min, x_max): mean_entry = current_mean dist_mean = 0 for c in range(5): @@ -117,7 +117,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, break # recompute means: means_list = [np.bincount(nearest_mean.ravel(), - image_yx[:, :, j].ravel()) for j in xrange(5)] + 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") diff --git a/skimage/segmentation/setup.py b/skimage/segmentation/setup.py index ec092ffe..b9e19078 100644 --- a/skimage/segmentation/setup.py +++ b/skimage/segmentation/setup.py @@ -11,8 +11,8 @@ def configuration(parent_package='', top_path=None): config = Configuration('segmentation', parent_package, top_path) - cython(['felzenszwalb_cy.pyx'], working_path=base_path) - config.add_extension('felzenszwalb_cy', sources=['felzenszwalb_cy.c'], + 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'],