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