diff --git a/skimage/color/__init__.py b/skimage/color/__init__.py
index 9202bee9..19620cb1 100644
--- a/skimage/color/__init__.py
+++ b/skimage/color/__init__.py
@@ -1,4 +1,5 @@
 from .colorconv import (convert_colorspace,
+                        guess_spatial_dimensions,
                         rgb2hsv,
                         hsv2rgb,
                         rgb2xyz,
@@ -45,6 +46,7 @@ from .colorlabel import color_dict, label2rgb
 
 
 __all__ = ['convert_colorspace',
+           'guess_spatial_dimensions',
            'rgb2hsv',
            'hsv2rgb',
            'rgb2xyz',
diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index 1ef5a6c9..d2b20316 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -49,6 +49,37 @@ from ..util import dtype
 from skimage._shared.utils import deprecated
 
 
+def guess_spatial_dimensions(image):
+    """Make an educated guess about whether an image has a channels dimension.
+
+    Parameters
+    ----------
+    image : ndarray
+        The input image.
+
+    Returns
+    -------
+    spatial_dims : int or None
+        The number of spatial dimensions of `image`. If ambiguous, the value
+        is `None`.
+
+    Raises
+    ------
+    ValueError
+        If the image array has less than two or more than four dimensions.
+    """
+    if image.ndim == 2:
+        return 2
+    if image.ndim == 3 and image.shape[-1] != 3:
+        return 3
+    if image.ndim == 3 and image.shape[-1] == 3:
+        return None
+    if image.ndim == 4 and image.shape[-1] == 3:
+        return 3
+    else:
+        raise ValueError("Expected 2D, 3D, or 4D array, got %iD." % image.ndim)
+
+
 @deprecated()
 def is_rgb(image):
     """Test whether the image is RGB or RGBA.
@@ -72,7 +103,7 @@ def is_gray(image):
         Input image.
 
     """
-    return np.squeeze(image).ndim == 2
+    return image.ndim in (2, 3) and not is_rgb(image)
 
 
 def convert_colorspace(arr, fromspace, tospace):
@@ -129,8 +160,9 @@ def _prepare_colorarray(arr):
     """
     arr = np.asanyarray(arr)
 
-    if arr.ndim != 3 or arr.shape[2] != 3:
-        msg = "the input array must be have a shape == (.,.,3))"
+    if arr.ndim not in [3, 4] or arr.shape[-1] != 3:
+        msg = ("the input array must be have a shape == (.., ..,[ ..,] 3)), " +
+               "got (" + (", ".join(map(str, arr.shape))) + ")")
         raise ValueError(msg)
 
     return dtype.img_as_float(arr)
@@ -413,12 +445,12 @@ def _convert(matrix, arr):
         The converted array.
     """
     arr = _prepare_colorarray(arr)
-    arr = np.swapaxes(arr, 0, 2)
+    arr = np.swapaxes(arr, 0, -1)
     oldshape = arr.shape
     arr = np.reshape(arr, (3, -1))
     out = np.dot(matrix, arr)
     out.shape = oldshape
-    out = np.swapaxes(out, 2, 0)
+    out = np.swapaxes(out, -1, 0)
 
     return np.ascontiguousarray(out)
 
@@ -473,17 +505,19 @@ def rgb2xyz(rgb):
     Parameters
     ----------
     rgb : array_like
-        The image in RGB format, in a 3-D array of shape (.., .., 3).
+        The image in RGB format, in a 3- or 4-D array of shape
+        (.., ..,[ ..,] 3).
 
     Returns
     -------
     out : ndarray
-        The image in XYZ format, in a 3-D array of shape (.., .., 3).
+        The image in XYZ format, in a 3- or 4-D array of shape
+        (.., ..,[ ..,] 3).
 
     Raises
     ------
     ValueError
-        If `rgb` is not a 3-D array of shape (.., .., 3).
+        If `rgb` is not a 3- or 4-D array of shape (.., ..,[ ..,] 3).
 
     Notes
     -----
@@ -628,23 +662,24 @@ def gray2rgb(image):
     Parameters
     ----------
     image : array_like
-        Input image of shape ``(M, N)``.
+        Input image of shape ``(M, N [, P])``.
 
     Returns
     -------
     rgb : ndarray
-        RGB image of shape ``(M, N, 3)``.
+        RGB image of shape ``(M, N, [, P], 3)``.
 
     Raises
     ------
     ValueError
-        If the input is not 2-dimensional.
+        If the input is not a 2- or 3-dimensional image.
 
     """
     if np.squeeze(image).ndim == 3 and image.shape[2] in (3, 4):
         return image
-    elif image.ndim == 2 or np.squeeze(image).ndim == 2:
-        return np.dstack((image, image, image))
+    elif is_gray(image):
+        image = image[..., np.newaxis]
+        return np.concatenate((image,)*3, axis=-1)
     else:
         raise ValueError("Input image expected to be RGB, RGBA or gray.")
 
@@ -655,17 +690,19 @@ def xyz2lab(xyz):
     Parameters
     ----------
     xyz : array_like
-        The image in XYZ format, in a 3-D array of shape (.., .., 3).
+        The image in XYZ format, in a 3- or 4-D array of shape
+        (.., ..,[ ..,] 3).
 
     Returns
     -------
     out : ndarray
-        The image in CIE-LAB format, in a 3-D array of shape (.., .., 3).
+        The image in CIE-LAB format, in a 3- or 4-D array of shape
+        (.., ..,[ ..,] 3).
 
     Raises
     ------
     ValueError
-        If `xyz` is not a 3-D array of shape (.., .., 3).
+        If `xyz` is not a 3-D array of shape (.., ..,[ ..,] 3).
 
     Notes
     -----
@@ -695,14 +732,14 @@ def xyz2lab(xyz):
     arr[mask] = np.power(arr[mask], 1. / 3.)
     arr[~mask] = 7.787 * arr[~mask] + 16. / 116.
 
-    x, y, z = arr[:, :, 0], arr[:, :, 1], arr[:, :, 2]
+    x, y, z = arr[..., 0], arr[..., 1], arr[..., 2]
 
     # Vector scaling
     L = (116. * y) - 16.
     a = 500.0 * (x - y)
     b = 200.0 * (y - z)
 
-    return np.dstack([L, a, b])
+    return np.concatenate([x[..., np.newaxis] for x in [L, a, b]], axis=-1)
 
 
 def lab2xyz(lab):
@@ -759,17 +796,19 @@ def rgb2lab(rgb):
     Parameters
     ----------
     rgb : array_like
-        The image in RGB format, in a 3-D array of shape (.., .., 3).
+        The image in RGB format, in a 3- or 4-D array of shape
+        (.., ..,[ ..,] 3).
 
     Returns
     -------
     out : ndarray
-        The image in Lab format, in a 3-D array of shape (.., .., 3).
+        The image in Lab format, in a 3- or 4-D array of shape
+        (.., ..,[ ..,] 3).
 
     Raises
     ------
     ValueError
-        If `rgb` is not a 3-D array of shape (.., .., 3).
+        If `rgb` is not a 3- or 4-D array of shape (.., ..,[ ..,] 3).
 
     Notes
     -----
diff --git a/skimage/color/tests/test_colorconv.py b/skimage/color/tests/test_colorconv.py
index 05ab6915..4fdaa4c8 100644
--- a/skimage/color/tests/test_colorconv.py
+++ b/skimage/color/tests/test_colorconv.py
@@ -33,7 +33,8 @@ from skimage.color import (rgb2hsv, hsv2rgb,
                            rgb2grey, gray2rgb,
                            xyz2lab, lab2xyz,
                            lab2rgb, rgb2lab,
-                           is_rgb, is_gray
+                           is_rgb, is_gray,
+                           guess_spatial_dimensions
                            )
 
 from skimage import data_dir, data
@@ -41,6 +42,19 @@ from skimage import data_dir, data
 import colorsys
 
 
+def test_guess_spatial_dimensions():
+    im1 = np.zeros((5, 5))
+    im2 = np.zeros((5, 5, 5))
+    im3 = np.zeros((5, 5, 3))
+    im4 = np.zeros((5, 5, 5, 3))
+    im5 = np.zeros((5,))
+    assert_equal(guess_spatial_dimensions(im1), 2)
+    assert_equal(guess_spatial_dimensions(im2), 3)
+    assert_equal(guess_spatial_dimensions(im3), None)
+    assert_equal(guess_spatial_dimensions(im4), 3)
+    assert_raises(ValueError, guess_spatial_dimensions, im5)
+
+
 class TestColorconv(TestCase):
 
     img_rgb = imread(os.path.join(data_dir, 'color.png'))
diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py
index c3aa1afc..107111bd 100644
--- a/skimage/segmentation/__init__.py
+++ b/skimage/segmentation/__init__.py
@@ -1,6 +1,6 @@
 from .random_walker_segmentation import random_walker
 from ._felzenszwalb import felzenszwalb
-from ._slic import slic
+from .slic_superpixels import slic
 from ._quickshift import quickshift
 from .boundaries import find_boundaries, visualize_boundaries, mark_boundaries
 from ._clear_border import clear_border
diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 9a5374d6..f6c6788c 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -2,140 +2,100 @@
 #cython: boundscheck=False
 #cython: nonecheck=False
 #cython: wraparound=False
+import collections as coll
 import numpy as np
 from time import time
 from scipy import ndimage
 
 cimport numpy as cnp
 
-from ..util import img_as_float
+from ..util import img_as_float, regular_grid
 from ..color import rgb2lab, gray2rgb
 
 
-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.
-
+def _slic_cython(double[:, :, :, ::1] image_zyx,
+                 long[:, :, ::1] nearest_mean,
+                 double[:, :, ::1] distance,
+                 double[:, ::1] means,
+                 float ratio, int max_iter, int n_segments):
+    """Helper function for SLIC segmentation.
+    
     Parameters
     ----------
-    image : (width, height [, 3]) ndarray
-        Input image.
-    n_segments : int, optional (default 100)
-        The (approximate) number of labels in the segmented output image.
-    ratio: float, optional (default 10)
-        Balances color-space proximity and image-space proximity.
-        Higher values give more weight to color-space.
-    max_iter : int, optional (default 10)
-        Maximum number of iterations of k-means.
-    sigma : float, optional (default 1)
-        Width of Gaussian smoothing kernel for preprocessing. Zero means no
-        smoothing.
-    convert2lab : bool, optional (default True)
-        Whether the input should be converted to Lab colorspace prior to
-        segmentation.  For this purpose, the input is assumed to be RGB. Highly
-        recommended.
+    image_zyx : 4D np.ndarray of double, shape (Z, Y, X, 6)
+        The image with embedded coordinates, that is, `image_zyx[i, j, k]` is
+        `array([i, j, k, r, g, b])` or `array([i, j, k, L, a, b])`, depending
+        on the colorspace.
+    nearest_mean : 3D np.ndarray of long, shape (Z, Y, X)
+        The (initially empty) label field.
+    distance : 3D np.ndarray of double, shape (Z, Y, X)
+        The (initially infinity) array of distances to the nearest centroid.
+    means : 2D np.ndarray of double, shape (n_segments, 6)
+        The centroids obtained by SLIC.
+    ratio : float
+        The ratio of xyz-space and colorspace in the clustering.
+    max_iter : int
+        The maximum number of k-means iterations.
+    n_segments : int
+        The approximate/desired number of segments.
 
     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.
-
-    Examples
-    --------
-    >>> from skimage.segmentation import slic
-    >>> from skimage.data import lena
-    >>> img = lena()
-    >>> segments = slic(img, n_segments=100, ratio=10)
-    >>> # Increasing the ratio parameter yields more square regions
-    >>> segments = slic(img, n_segments=100, ratio=20)
+    nearest_mean : 3D np.ndarray of long, shape (Z, Y, X)
+        The label field/superpixels found by SLIC.
     """
-    if image.ndim == 2:
-        image = gray2rgb(image)
-    if image.ndim != 3 or image.shape[2] != 3:
-        ValueError("Only 1- or 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 Py_ssize_t height, width
-    height, width = image.shape[:2]
+    cdef Py_ssize_t depth, height, width
+    depth, height, width = (image_zyx.shape[0], image_zyx.shape[1],
+                            image_zyx.shape[2])
     # approximate grid size for desired n_segments
-    cdef Py_ssize_t step = int(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]
+    cdef Py_ssize_t step_z, step_y, step_x
+    slices = regular_grid((depth, height, width), n_segments)
+    step_z, step_y, step_x = [int(s.step) for s in slices]
 
-    means_color = np.zeros((means_y.shape[0], means_y.shape[1], 3))
-    cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] means \
-            = np.dstack([means_y, means_x, means_color]).reshape(-1, 5)
-    cdef cnp.float_t* current_mean
-    cdef cnp.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 cnp.ndarray[dtype=cnp.float_t, ndim=3] image_yx \
-            = np.dstack([grid_y, grid_x, image / ratio]).copy("C")
-    cdef Py_ssize_t i, k, x, y, x_min, x_max, y_min, y_max, changes
+    cdef Py_ssize_t i, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max, \
+            changes
     cdef double dist_mean
 
-    cdef cnp.ndarray[dtype=cnp.intp_t, ndim=2] nearest_mean \
-            = np.zeros((height, width), dtype=np.intp)
-    cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] distance \
-            = np.empty((height, width))
-    cdef cnp.float_t* image_p = <cnp.float_t*> image_yx.data
-    cdef cnp.float_t* distance_p = <cnp.float_t*> distance.data
-    cdef cnp.float_t* current_distance
-    cdef cnp.float_t* current_pixel
     cdef double tmp
     for i in range(max_iter):
-        distance.fill(np.inf)
         changes = 0
-        current_mean = <cnp.float_t*> means.data
+        distance[:, :, :] = np.inf
         # 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
+            z_min = int(max(means[k, 0] - 2 * step_z, 0))
+            z_max = int(min(means[k, 0] + 2 * step_z, depth))
+            y_min = int(max(means[k, 1] - 2 * step_y, 0))
+            y_max = int(min(means[k, 1] + 2 * step_y, height))
+            x_min = int(max(means[k, 2] - 2 * step_x, 0))
+            x_max = int(min(means[k, 2] + 2 * step_x, width))
+            for z in range(z_min, z_max):
+                for y in range(y_min, y_max):
+                    for x in range(x_min, x_max):
+                        dist_mean = 0
+                        for c in range(6):
+                            # you would think the compiler can optimize the
+                            # squaring itself. mine can't (with O2)
+                            tmp = image_zyx[z, y, x, c] - means[k, c]
+                            dist_mean += tmp * tmp
+                        # some precision issue here. Doesnt work if testing ">"
+                        if distance[z, y, x] - dist_mean > 1e-10:
+                            nearest_mean[z, y, x] = k
+                            distance[z, y, x] = dist_mean
+                            changes = 1
         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())
+        nearest_mean_ravel = np.asarray(nearest_mean).ravel()
+        means_list = []
+        for j in range(6):
+            image_zyx_ravel = np.ascontiguousarray(image_zyx[:, :, :, j]).ravel()
+            means_list.append(np.bincount(nearest_mean_ravel,
+                                          image_zyx_ravel))
+        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
+    return np.ascontiguousarray(nearest_mean)
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
new file mode 100644
index 00000000..2caf8ade
--- /dev/null
+++ b/skimage/segmentation/slic_superpixels.py
@@ -0,0 +1,136 @@
+# coding=utf-8
+
+import collections as coll
+import numpy as np
+from scipy import ndimage
+import warnings
+
+from ..util import img_as_float, regular_grid
+from ..color import rgb2lab, gray2rgb, guess_spatial_dimensions
+from ._slic import _slic_cython
+
+
+def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
+         multichannel=None, convert2lab=True):
+    """Segments image using k-means clustering in Color-(x,y) space.
+
+    Parameters
+    ----------
+    image : (width, height [, depth] [, 3]) ndarray
+        Input image, which can be 2D or 3D, and grayscale or multi-channel
+        (see `multichannel` parameter).
+    n_segments : int, optional (default: 100)
+        The (approximate) number of labels in the segmented output image.
+    ratio: float, optional (default: 10)
+        Balances color-space proximity and image-space proximity.
+        Higher values give more weight to color-space.
+    max_iter : int, optional (default: 10)
+        Maximum number of iterations of k-means.
+    sigma : float, optional (default: 1)
+        Width of Gaussian smoothing kernel for preprocessing. Zero means no
+        smoothing.
+    multichannel : bool, optional (default: None)
+        Whether the last axis of the image is to be interpreted as multiple
+        channels. Only 3 channels are supported. If `None`, the function will
+        attempt to guess this, and raise a warning if ambiguous, when the
+        array has shape (M, N, 3).
+    convert2lab : bool, optional (default: True)
+        Whether the input should be converted to Lab colorspace prior to
+        segmentation.  For this purpose, the input is assumed to be RGB. Highly
+        recommended.
+
+    Returns
+    -------
+    segment_mask : (width, height) ndarray
+        Integer mask indicating segment labels.
+
+    Raises
+    ------
+    ValueError
+        If:
+            - the image dimension is not 2 or 3 and `multichannel == False`, OR
+            - the image dimension is not 3 or 4 and `multichannel == True`, OR
+            - `multichannel == True` and the length of the last dimension of
+            the image is not 3, OR
+
+    Notes
+    -----
+    If `sigma > 0` as is default, the image is smoothed using a Gaussian kernel
+    prior to segmentation.
+
+    The image is rescaled to be in [0, 1] prior to processing.
+
+    Images of shape (M, N, 3) are interpreted as 2D RGB images by default. To
+    interpret them as 3D with the last dimension having length 3, use
+    `multichannel=False`.
+
+    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.
+
+    Examples
+    --------
+    >>> from skimage.segmentation import slic
+    >>> from skimage.data import lena
+    >>> img = lena()
+    >>> segments = slic(img, n_segments=100, ratio=10)
+    >>> # Increasing the ratio parameter yields more square regions
+    >>> segments = slic(img, n_segments=100, ratio=20)
+    """
+    spatial_dims = guess_spatial_dimensions(image)
+    if spatial_dims is None and multichannel is None:
+        msg = ("Images with dimensions (M, N, 3) are interpreted as 2D+RGB" +
+                   " by default. Use `multichannel=False` to interpret as " +
+                   " 3D image with last dimension of length 3.")
+        warnings.warn(RuntimeWarning(msg))
+        multichannel = True
+    elif multichannel is None:
+        multichannel = (spatial_dims + 1 == image.ndim)
+    if ((not multichannel and image.ndim not in [2, 3]) or
+            (multichannel and image.ndim not in [3, 4]) or
+            (multichannel and image.shape[-1] != 3)):
+        ValueError("Only 1- or 3-channel 2- or 3-D images are supported.")
+    image = img_as_float(image)
+    if not multichannel:
+        image = gray2rgb(image)
+    if image.ndim == 3:
+        # See 2D RGB image as 3D RGB image with Z = 1
+        image = image[np.newaxis, ...]
+    if not isinstance(sigma, coll.Iterable):
+        sigma = np.array([sigma, sigma, sigma, 0])
+    if (sigma > 0).any():
+        image = ndimage.gaussian_filter(image, sigma)
+    if convert2lab:
+        image = rgb2lab(image)
+
+    # initialize on grid:
+    depth, height, width = image.shape[:3]
+    # approximate grid size for desired n_segments
+    grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
+    slices = regular_grid(image.shape[:3], n_segments)
+    step_z, step_y, step_x = [int(s.step) for s in slices]
+    means_z = grid_z[slices]
+    means_y = grid_y[slices]
+    means_x = grid_x[slices]
+
+    means_color = np.zeros(means_z.shape + (3,))
+    means = np.concatenate([means_z[..., np.newaxis], means_y[..., np.newaxis],
+                            means_x[..., np.newaxis], means_color
+                           ], axis=-1).reshape(-1, 6)
+    means = np.ascontiguousarray(means)
+    # 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(max((step_z, step_y, step_x)))) ** 2
+    image_zyx = np.concatenate([grid_z[..., np.newaxis],
+                                grid_y[..., np.newaxis],
+                                grid_x[..., np.newaxis],
+                                image / ratio], axis=-1).copy("C")
+    nearest_mean = np.zeros((depth, height, width), dtype=np.intp)
+    distance = np.empty((depth, height, width), dtype=np.float)
+    segment_map = _slic_cython(image_zyx, nearest_mean, distance, means,
+                               ratio, max_iter, n_segments)
+    if segment_map.shape[0] == 1:
+        segment_map = segment_map[0]
+    return segment_map
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 89dee59b..d0539cd2 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -1,9 +1,11 @@
+import itertools as it
+import warnings
 import numpy as np
 from numpy.testing import assert_equal, assert_array_equal
 from skimage.segmentation import slic
 
 
-def test_color():
+def test_color_2d():
     rnd = np.random.RandomState(0)
     img = np.zeros((20, 21, 3))
     img[:10, :10, 0] = 1
@@ -12,7 +14,9 @@ def test_color():
     img += 0.01 * rnd.normal(size=img.shape)
     img[img > 1] = 1
     img[img < 0] = 0
-    seg = slic(img, sigma=0, n_segments=4)
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        seg = slic(img, sigma=0, n_segments=4)
 
     # we expect 4 segments
     assert_equal(len(np.unique(seg)), 4)
@@ -21,7 +25,8 @@ def test_color():
     assert_array_equal(seg[:10, 10:], 1)
     assert_array_equal(seg[10:, 10:], 3)
 
-def test_gray():
+
+def test_gray_2d():
     rnd = np.random.RandomState(0)
     img = np.zeros((20, 21))
     img[:10, :10] = 0.33
@@ -30,7 +35,7 @@ def test_gray():
     img += 0.0033 * rnd.normal(size=img.shape)
     img[img > 1] = 1
     img[img < 0] = 0
-    seg = slic(img, sigma=0, n_segments=4, ratio=50.0)
+    seg = slic(img, sigma=0, n_segments=4, ratio=20.0, multichannel=False)
 
     assert_equal(len(np.unique(seg)), 4)
     assert_array_equal(seg[:10, :10], 0)
@@ -38,6 +43,49 @@ def test_gray():
     assert_array_equal(seg[:10, 10:], 1)
     assert_array_equal(seg[10:, 10:], 3)
 
+
+def test_color_3d():
+    rnd = np.random.RandomState(0)
+    img = np.zeros((20, 21, 22, 3))
+    slices = []
+    for dim_size in img.shape[:-1]:
+        midpoint = dim_size // 2
+        slices.append((slice(None, midpoint), slice(midpoint, None)))
+    slices = list(it.product(*slices))
+    colors = list(it.product(*(([0, 1],) * 3)))
+    for s, c in zip(slices, colors):
+        img[s] = c
+    img += 0.01 * rnd.normal(size=img.shape)
+    img[img > 1] = 1
+    img[img < 0] = 0
+    seg = slic(img, sigma=0, n_segments=8)
+
+    assert_equal(len(np.unique(seg)), 8)
+    for s, c in zip(slices, range(8)):
+        assert_array_equal(seg[s], c)
+
+
+def test_gray_3d():
+    rnd = np.random.RandomState(0)
+    img = np.zeros((20, 21, 22))
+    slices = []
+    for dim_size in img.shape:
+        midpoint = dim_size // 2
+        slices.append((slice(None, midpoint), slice(midpoint, None)))
+    slices = list(it.product(*slices))
+    shades = np.arange(0, 1.000001, 1.0/7)
+    for s, sh in zip(slices, shades):
+        img[s] = sh
+    img += 0.001 * rnd.normal(size=img.shape)
+    img[img > 1] = 1
+    img[img < 0] = 0
+    seg = slic(img, sigma=0, n_segments=8, ratio=20.0, multichannel=False)
+
+    assert_equal(len(np.unique(seg)), 8)
+    for s, c in zip(slices, range(8)):
+        assert_array_equal(seg[s], c)
+
+
 if __name__ == '__main__':
     from numpy import testing
     testing.run_module_suite()
diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index a4274484..7afcf54a 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -11,6 +11,7 @@ if chk < 18:    # Use internal version for numpy versions < 1.8.x
 else:
     from numpy import pad
 del numpy, ver, chk
+from ._regular_grid import regular_grid
 
 
 __all__ = ['img_as_float',
@@ -22,4 +23,5 @@ __all__ = ['img_as_float',
            'view_as_blocks',
            'view_as_windows',
            'pad',
-           'random_noise']
+           'random_noise',
+           'regular_grid']
diff --git a/skimage/util/_regular_grid.py b/skimage/util/_regular_grid.py
new file mode 100644
index 00000000..e304be20
--- /dev/null
+++ b/skimage/util/_regular_grid.py
@@ -0,0 +1,72 @@
+import numpy as np
+
+
+def regular_grid(ar_shape, n_points):
+    """Find `n_points` regularly spaced along `ar_shape`.
+    
+    The returned points (as slices) should be as close to cubically-spaced as
+    possible. Essentially, the points are spaced by the Nth root of the input
+    array size, where N is the number of dimensions. However, if an array
+    dimension cannot fit a full step size, it is "discarded", and the
+    computation is done for only the remaining dimensions.
+
+    Parameters
+    ----------
+    ar_shape : array-like of ints
+        The shape of the space embedding the grid. `len(ar_shape)` is the
+        number of dimensions.
+    n_points : int
+        The (approximate) number of points to embed in the space.
+
+    Returns
+    -------
+    slices : list of slice objects
+        A slice along each dimension of `ar_shape`, such that the intersection
+        of all the slices give the coordinates of regularly spaced points.
+
+    Examples
+    --------
+    >>> ar = np.zeros((20, 40))
+    >>> g = regular_grid(ar.shape, 8)
+    >>> g
+    [slice(5.0, None, 10.0), slice(5.0, None, 10.0)]
+    >>> ar[g] = 1
+    >>> ar.sum()
+    8.0
+    >>> ar = np.zeros((20, 40))
+    >>> g = regular_grid(ar.shape, 32)
+    >>> g
+    [slice(2.0, None, 5.0), slice(2.0, None, 5.0)]
+    >>> ar[g] = 1
+    >>> ar.sum()
+    32.0
+    >>> ar = np.zeros((3, 20, 40))
+    >>> g = regular_grid(ar.shape, 8)
+    >>> g
+    [slice(1.0, None, 3.0), slice(5.0, None, 10.0), slice(5.0, None, 10.0)]
+    >>> ar[g] = 1
+    >>> ar.sum()
+    8.0
+    """
+    ar_shape = np.asanyarray(ar_shape)
+    ndim = len(ar_shape)
+    unsort_dim_idxs = np.argsort(np.argsort(ar_shape))
+    sorted_dims = np.sort(ar_shape)
+    space_size = float(np.prod(ar_shape))
+    if space_size <= n_points:
+        return [slice(None)] * ndim
+    stepsizes = (space_size / n_points) ** (1.0 / ndim) * np.ones(ndim)
+    if (sorted_dims < stepsizes).any():
+        for dim in range(ndim):
+            stepsizes[dim] = sorted_dims[dim]
+            space_size = float(np.prod(sorted_dims[dim+1:]))
+            stepsizes[dim+1:] = ((space_size / n_points) **
+                                 (1.0 / (ndim - dim - 1)))
+            if (sorted_dims >= stepsizes).all():
+                break
+    starts = stepsizes // 2
+    stepsizes = np.round(stepsizes)
+    slices = [slice(start, None, step) for 
+              start, step in zip(starts, stepsizes)]
+    slices = [slices[i] for i in unsort_dim_idxs]
+    return slices
diff --git a/skimage/util/tests/test_regular_grid.py b/skimage/util/tests/test_regular_grid.py
new file mode 100644
index 00000000..61736a76
--- /dev/null
+++ b/skimage/util/tests/test_regular_grid.py
@@ -0,0 +1,40 @@
+import numpy as np
+from numpy.testing import assert_equal
+from skimage.util import regular_grid
+
+
+def test_regular_grid_full():
+    ar = np.zeros((2, 2))
+    g = regular_grid(ar, 25)
+    assert_equal(g, [slice(None, None, None), slice(None, None, None)])
+    ar[g] = 1
+    assert_equal(ar.size, ar.sum())
+
+
+def test_regular_grid_2d_8():
+    ar = np.zeros((20, 40))
+    g = regular_grid(ar.shape, 8)
+    assert_equal(g, [slice(5.0, None, 10.0), slice(5.0, None, 10.0)])
+    ar[g] = 1
+    assert_equal(ar.sum(), 8)
+
+
+def test_regular_grid_2d_32():
+    ar = np.zeros((20, 40))
+    g = regular_grid(ar.shape, 32)
+    assert_equal(g, [slice(2.0, None, 5.0), slice(2.0, None, 5.0)])
+    ar[g] = 1
+    assert_equal(ar.sum(), 32)
+
+
+def test_regular_grid_3d_8():
+    ar = np.zeros((3, 20, 40))
+    g = regular_grid(ar.shape, 8)
+    assert_equal(g, [slice(1.0, None, 3.0), slice(5.0, None, 10.0),
+                     slice(5.0, None, 10.0)])
+    ar[g] = 1
+    assert_equal(ar.sum(), 8)
+
+
+if __name__ == '__main__':
+    np.testing.run_module_suite()