From 19f264bc8b7092fc6e3b8f86181fbf7679b24d52 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 11 Mar 2013 21:03:12 +1100
Subject: [PATCH 01/39] Make color functions used by SLIC 3D-aware

---
 skimage/color/colorconv.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index 1ef5a6c9..6789dbbc 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -59,7 +59,7 @@ def is_rgb(image):
         Input image.
 
     """
-    return (image.ndim == 3 and image.shape[2] in (3, 4))
+    return (image.ndim in (3, 4) and image.shape[-1] in (3, 4))
 
 
 @deprecated()
@@ -72,7 +72,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):
@@ -628,23 +628,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.")
 

From 0132998ca111b4cd35620d05d27a82877636c02c Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 12 Mar 2013 00:28:49 +1100
Subject: [PATCH 02/39] Make more colorconv functions 3D aware

---
 skimage/color/colorconv.py | 36 +++++++++++++++++++++---------------
 1 file changed, 21 insertions(+), 15 deletions(-)

diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index 6789dbbc..be2f1519 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -129,8 +129,8 @@ 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))"
         raise ValueError(msg)
 
     return dtype.img_as_float(arr)
@@ -413,12 +413,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 +473,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
     -----
@@ -656,17 +658,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
     -----
@@ -696,14 +700,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(map(lambda x: x[..., np.newaxis], [L, a, b]), -1)
 
 
 def lab2xyz(lab):
@@ -760,17 +764,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
     -----

From a2e32cc90c5271308ec32ee588cfeff5ae6fd445 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Thu, 14 Mar 2013 01:56:14 +1100
Subject: [PATCH 03/39] Add initial 3D modifications (not working)

---
 skimage/segmentation/_slic.pyx          | 152 +++++++++++++++---------
 skimage/segmentation/tests/test_slic.py |   2 +-
 2 files changed, 98 insertions(+), 56 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 9a5374d6..35a55804 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -2,6 +2,7 @@
 #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
@@ -13,24 +14,28 @@ from ..color import rgb2lab, gray2rgb
 
 
 def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
-         convert2lab=True):
+         multichannel=True, convert2lab=True):
     """Segments image using k-means clustering in Color-(x,y) space.
 
     Parameters
     ----------
-    image : (width, height [, 3]) ndarray
-        Input image.
-    n_segments : int, optional (default 100)
+    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)
+    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)
+    max_iter : int, optional (default: 10)
         Maximum number of iterations of k-means.
-    sigma : float, optional (default 1)
+    sigma : float, optional (default: 1)
         Width of Gaussian smoothing kernel for preprocessing. Zero means no
         smoothing.
-    convert2lab : bool, optional (default True)
+    multichannel : bool, optional (default: True)
+        Whether the last axis of the image is to be interpreted as multiple
+        channels. Only 3 channels are supported.
+    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.
@@ -40,9 +45,19 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     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.
+
     Notes
     -----
-    The image is smoothed using a Gaussian kernel prior to segmentation.
+    The image is optionally smoothed using a Gaussian kernel prior to
+    segmentation.
 
     References
     ----------
@@ -59,42 +74,64 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     >>> # Increasing the ratio parameter yields more square regions
     >>> segments = slic(img, n_segments=100, ratio=20)
     """
-    if image.ndim == 2:
+    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.")
+    if image.ndim in [2, 3] and not multichannel:
         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 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(img_as_float(image), sigma)
     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.shape[:3]
     # 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 = int(np.ceil(
+                              (depth * height * width / n_segments) **
+                              (1.0/3)))
+    grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
+    means_z = grid_z[::step, ::step, ::step]
+    means_y = grid_y[::step, ::step, ::step]
+    means_x = grid_x[::step, ::step, ::step]
 
-    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)
+    means_color = np.zeros(means_z.shape + (3,))
+    cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] means = \
+            np.concatenate([
+                            means_z[..., np.newaxis],
+                            means_y[..., np.newaxis],
+                            means_x[..., np.newaxis],
+                            means_color
+                           ], axis=-1).reshape(-1, 6)
     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 cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx \
+            = np.concatenate([
+                              grid_y[..., np.newaxis],
+                              grid_x[..., np.newaxis],
+                              grid_z[..., np.newaxis],
+                              image / ratio
+                             ], axis=-1).copy("C")
+    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.ndarray[dtype=cnp.intp_t, ndim=3] nearest_mean \
+            = np.zeros((depth, height, width), dtype=np.intp)
+    cdef cnp.ndarray[dtype=cnp.float_t, ndim=3] distance \
+            = np.empty((depth, height, width))
+    cdef cnp.float_t* image_p = <cnp.float_t*> image_zyx.data
     cdef cnp.float_t* distance_p = <cnp.float_t*> distance.data
     cdef cnp.float_t* current_distance
     cdef cnp.float_t* current_pixel
@@ -106,35 +143,40 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
         # 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(current_mean[0] - 2 * step, 0))
+            z_max = int(min(current_mean[0] + 2 * step, depth))
+            y_min = int(max(current_mean[1] - 2 * step, 0))
+            y_max = int(min(current_mean[1] + 2 * step, height))
+            x_min = int(max(current_mean[2] - 2 * step, 0))
+            x_max = int(min(current_mean[2] + 2 * step, width))
+            for z in range(z_min, z_max):
+                for y in range(y_min, y_max):
+                    current_pixel = \
+                            &image_p[5 * ((z * height + y) * width + x_min)]
+                    current_distance = \
+                            &distance_p[(z * height + 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[z, y, x] = k
+                            current_distance[0] = dist_mean
+                            changes += 1
+                        current_distance += 1
+            current_mean += 6
         if changes == 0:
             break
         # recompute means:
         means_list = [np.bincount(nearest_mean.ravel(),
-                      image_yx[:, :, j].ravel()) for j in range(5)]
+                      image_zyx[:, :, :, j].ravel()) for j in range(6)]
         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/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 89dee59b..b080e378 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -30,7 +30,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=50.0, multichannel=False)
 
     assert_equal(len(np.unique(seg)), 4)
     assert_array_equal(seg[:10, :10], 0)

From bb8cfea8c63822565987a623b8d66e7a405a8e6c Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Thu, 14 Mar 2013 18:13:49 +1100
Subject: [PATCH 04/39] Add function to calculate regularly-spaced grid in nD

---
 skimage/util/__init__.py     |  2 ++
 skimage/util/regular_grid.py | 44 ++++++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+)
 create mode 100644 skimage/util/regular_grid.py

diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index a4274484..474b2966 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',
@@ -23,3 +24,4 @@ __all__ = ['img_as_float',
            'view_as_windows',
            'pad',
            'random_noise']
+           'regular_grid']
diff --git a/skimage/util/regular_grid.py b/skimage/util/regular_grid.py
new file mode 100644
index 00000000..40980491
--- /dev/null
+++ b/skimage/util/regular_grid.py
@@ -0,0 +1,44 @@
+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.
+
+    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.
+    """
+    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 = np.floor(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

From 620210025f5dd2d363eb6815776220ee56626d3e Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Thu, 14 Mar 2013 18:33:01 +1100
Subject: [PATCH 05/39] Modify SLIC to allow uneven step sizes

---
 skimage/segmentation/_slic.pyx | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 35a55804..3ffbb5f6 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -9,7 +9,7 @@ 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
 
 
@@ -94,13 +94,13 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     cdef Py_ssize_t depth, height, width
     depth, height, width = image.shape[:3]
     # approximate grid size for desired n_segments
-    cdef Py_ssize_t step = int(np.ceil(
-                              (depth * height * width / n_segments) **
-                              (1.0/3)))
+    cdef Py_ssize_t step_z, step_y, step_x
     grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
-    means_z = grid_z[::step, ::step, ::step]
-    means_y = grid_y[::step, ::step, ::step]
-    means_x = grid_x[::step, ::step, ::step]
+    slices = regular_grid(image.shape, 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,))
     cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] means = \
@@ -115,7 +115,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     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
+    ratio = (ratio / float(max((step_z, step_y, step_x)))) ** 2
     cdef cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx \
             = np.concatenate([
                               grid_y[..., np.newaxis],
@@ -143,12 +143,12 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
         # assign pixels to means
         for k in range(n_means):
             # compute windows:
-            z_min = int(max(current_mean[0] - 2 * step, 0))
-            z_max = int(min(current_mean[0] + 2 * step, depth))
-            y_min = int(max(current_mean[1] - 2 * step, 0))
-            y_max = int(min(current_mean[1] + 2 * step, height))
-            x_min = int(max(current_mean[2] - 2 * step, 0))
-            x_max = int(min(current_mean[2] + 2 * step, width))
+            z_min = int(max(current_mean[0] - 2 * step_z, 0))
+            z_max = int(min(current_mean[0] + 2 * step_z, depth))
+            y_min = int(max(current_mean[1] - 2 * step_y, 0))
+            y_max = int(min(current_mean[1] + 2 * step_y, height))
+            x_min = int(max(current_mean[2] - 2 * step_x, 0))
+            x_max = int(min(current_mean[2] + 2 * step_x, width))
             for z in range(z_min, z_max):
                 for y in range(y_min, y_max):
                     current_pixel = \

From 30636428db4e0749b612645cb0e46bfc6a9a3373 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Thu, 14 Mar 2013 18:52:30 +1100
Subject: [PATCH 06/39] Add more descriptive error message to
 _prepare_colorarray

---
 skimage/color/colorconv.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index be2f1519..540c5a75 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -130,7 +130,8 @@ def _prepare_colorarray(arr):
     arr = np.asanyarray(arr)
 
     if arr.ndim not in [3, 4] or arr.shape[-1] != 3:
-        msg = "the input array must be have a shape == (.., ..,[ ..,] 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)

From 6dc8e6300bb8e0c7aa11b67575605e5e1690912b Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Wed, 3 Apr 2013 17:39:56 +1100
Subject: [PATCH 07/39] Separate inner loop of SLIC computation

---
 skimage/segmentation/__init__.py        |   2 +-
 skimage/segmentation/_slic.pyx          | 115 +++---------------------
 skimage/segmentation/slic.py            | 114 +++++++++++++++++++++++
 skimage/segmentation/tests/test_slic.py |   4 +-
 4 files changed, 127 insertions(+), 108 deletions(-)
 create mode 100644 skimage/segmentation/slic.py

diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py
index c3aa1afc..ae1bf074 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 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 3ffbb5f6..20e0001a 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -13,124 +13,29 @@ 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,
-         multichannel=True, 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: True)
-        Whether the last axis of the image is to be interpreted as multiple
-        channels. Only 3 channels are supported.
-    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.
-
-    Notes
-    -----
-    The image is optionally 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)
-    """
-    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.")
-    if image.ndim in [2, 3] and 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(img_as_float(image), sigma)
-    if convert2lab:
-        image = rgb2lab(image)
+def _slic_cython(cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx,
+                 cnp.ndarray[dtype=cnp.intp_t, ndim=3] nearest_mean,
+                 cnp.ndarray[dtype=cnp.float_t, ndim=3] distance,
+                 cnp.ndarray[dtype=cnp.float_t, ndim=2] means,
+                 float ratio, int max_iter, int n_segments):
+    """Helper function for SLIC segmentation."""
 
     # initialize on grid:
     cdef Py_ssize_t depth, height, width
-    depth, height, width = image.shape[:3]
+    depth, height, width = image_zyx.shape[0], image_zyx[1], image_zyx[2]
     # approximate grid size for desired n_segments
     cdef Py_ssize_t step_z, step_y, step_x
     grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
-    slices = regular_grid(image.shape, n_segments)
+    slices = regular_grid((depth, height, width), 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,))
-    cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] means = \
-            np.concatenate([
-                            means_z[..., np.newaxis],
-                            means_y[..., np.newaxis],
-                            means_x[..., np.newaxis],
-                            means_color
-                           ], axis=-1).reshape(-1, 6)
     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(max((step_z, step_y, step_x)))) ** 2
-    cdef cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx \
-            = np.concatenate([
-                              grid_y[..., np.newaxis],
-                              grid_x[..., np.newaxis],
-                              grid_z[..., np.newaxis],
-                              image / ratio
-                             ], axis=-1).copy("C")
     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=3] nearest_mean \
-            = np.zeros((depth, height, width), dtype=np.intp)
-    cdef cnp.ndarray[dtype=cnp.float_t, ndim=3] distance \
-            = np.empty((depth, height, width))
     cdef cnp.float_t* image_p = <cnp.float_t*> image_zyx.data
     cdef cnp.float_t* distance_p = <cnp.float_t*> distance.data
     cdef cnp.float_t* current_distance
@@ -152,13 +57,13 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
             for z in range(z_min, z_max):
                 for y in range(y_min, y_max):
                     current_pixel = \
-                            &image_p[5 * ((z * height + y) * width + x_min)]
+                            &image_p[6 * ((z * height + y) * width + x_min)]
                     current_distance = \
                             &distance_p[(z * height + y) * width + x_min]
                     for x in range(x_min, x_max):
                         mean_entry = current_mean
                         dist_mean = 0
-                        for c in range(5):
+                        for c in range(6):
                             # you would think the compiler can optimize the
                             # squaring itself. mine can't (with O2)
                             tmp = current_pixel[0] - mean_entry[0]
diff --git a/skimage/segmentation/slic.py b/skimage/segmentation/slic.py
new file mode 100644
index 00000000..b5d28f11
--- /dev/null
+++ b/skimage/segmentation/slic.py
@@ -0,0 +1,114 @@
+import collections as coll
+import numpy as np
+from scipy import ndimage
+
+from ..util import img_as_float, regular_grid
+from ..color import rgb2lab, gray2rgb
+from ._slic import _slic_cython
+
+
+def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
+         multichannel=True, 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: True)
+        Whether the last axis of the image is to be interpreted as multiple
+        channels. Only 3 channels are supported.
+    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.
+
+    Notes
+    -----
+    The image is optionally 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)
+    """
+    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.")
+    if image.ndim in [2, 3] and 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(img_as_float(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)
+    # 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_y[..., np.newaxis],
+                                grid_x[..., np.newaxis],
+                                grid_z[..., 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 b080e378..58fda099 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -12,7 +12,7 @@ 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)
+    seg = slic(img, sigma=0, n_segments=4)[0]
 
     # we expect 4 segments
     assert_equal(len(np.unique(seg)), 4)
@@ -30,7 +30,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, multichannel=False)
+    seg = slic(img, sigma=0, n_segments=4, ratio=50.0, multichannel=False)[0]
 
     assert_equal(len(np.unique(seg)), 4)
     assert_array_equal(seg[:10, :10], 0)

From 8e4ab32ed9c5ba48802aadb653ab4055c1826345 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 13 May 2013 17:07:42 +1000
Subject: [PATCH 08/39] Bug fixes: concatenation order and shape assignments

---
 skimage/segmentation/_slic.pyx | 3 ++-
 skimage/segmentation/slic.py   | 4 ++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 20e0001a..d668f084 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -22,7 +22,8 @@ def _slic_cython(cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx,
 
     # initialize on grid:
     cdef Py_ssize_t depth, height, width
-    depth, height, width = image_zyx.shape[0], image_zyx[1], image_zyx[2]
+    shape = image_zyx.shape
+    depth, height, width = shape[0], shape[1], shape[2]
     # approximate grid size for desired n_segments
     cdef Py_ssize_t step_z, step_y, step_x
     grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
diff --git a/skimage/segmentation/slic.py b/skimage/segmentation/slic.py
index b5d28f11..9cf95f11 100644
--- a/skimage/segmentation/slic.py
+++ b/skimage/segmentation/slic.py
@@ -101,9 +101,9 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     # 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_y[..., np.newaxis],
+    image_zyx = np.concatenate([grid_z[..., np.newaxis],
+                                grid_y[..., np.newaxis],
                                 grid_x[..., np.newaxis],
-                                grid_z[..., 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)

From e9aa78b937b1d6dd18c895b346b2e6eb4a23cee0 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 14 May 2013 03:06:01 +1000
Subject: [PATCH 09/39] Bug fix: don't add singleton dimension to 3D gray
 images

---
 skimage/segmentation/slic.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/slic.py b/skimage/segmentation/slic.py
index 9cf95f11..d039d5ed 100644
--- a/skimage/segmentation/slic.py
+++ b/skimage/segmentation/slic.py
@@ -3,7 +3,7 @@ import numpy as np
 from scipy import ndimage
 
 from ..util import img_as_float, regular_grid
-from ..color import rgb2lab, gray2rgb
+from ..color import rgb2lab, gray2rgb, is_rgb
 from ._slic import _slic_cython
 
 
@@ -72,9 +72,9 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
             (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.")
-    if image.ndim in [2, 3] and not multichannel:
+    if not multichannel:
         image = gray2rgb(image)
-    if image.ndim == 3:
+    if image.ndim == 3 and is_rgb(image):
         # See 2D RGB image as 3D RGB image with Z = 1
         image = image[np.newaxis, ...]
     if not isinstance(sigma, coll.Iterable):

From 57cc86d7c81565d8f51037ebc178e90bc8981052 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 14 May 2013 03:07:09 +1000
Subject: [PATCH 10/39] Bug fix: remove unnecessary __get__ in test_slic

---
 skimage/segmentation/tests/test_slic.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 58fda099..1f5f0c42 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -12,7 +12,7 @@ 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)[0]
+    seg = slic(img, sigma=0, n_segments=4)
 
     # we expect 4 segments
     assert_equal(len(np.unique(seg)), 4)
@@ -30,7 +30,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, multichannel=False)[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)

From 215439b43c27271c95fc208bf683a19619c81b8d Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 14 May 2013 03:07:28 +1000
Subject: [PATCH 11/39] Add 3D slic tests (gray not working yet)

---
 skimage/segmentation/tests/test_slic.py | 49 ++++++++++++++++++++++++-
 1 file changed, 47 insertions(+), 2 deletions(-)

diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 1f5f0c42..9a277970 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -1,9 +1,10 @@
+import itertools as it
 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
@@ -21,7 +22,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
@@ -38,6 +40,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[:-1]:
+        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=40.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()

From 5c4d0218cea97cfbe761906e7fda38b2b7381330 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 14 May 2013 03:24:27 +1000
Subject: [PATCH 12/39] bug fix: 3d test working

---
 skimage/segmentation/tests/test_slic.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 9a277970..1018969a 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -66,7 +66,7 @@ def test_gray_3d():
     rnd = np.random.RandomState(0)
     img = np.zeros((20, 21, 22))
     slices = []
-    for dim_size in img.shape[:-1]:
+    for dim_size in img.shape:
         midpoint = dim_size // 2
         slices.append((slice(None, midpoint), slice(midpoint, None)))
     slices = list(it.product(*slices))
@@ -76,7 +76,7 @@ def test_gray_3d():
     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=40.0, multichannel=False)
+    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)):

From 332ab0449be9bff38c2edece247a1f09ef5a28de Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 14 May 2013 13:15:29 +1000
Subject: [PATCH 13/39] Improve PEP8 and Python 3 compliance

Thanks to @JDWarner and @ahojnnes for the input.
---
 skimage/util/regular_grid.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/util/regular_grid.py b/skimage/util/regular_grid.py
index 40980491..ab1b12dc 100644
--- a/skimage/util/regular_grid.py
+++ b/skimage/util/regular_grid.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def regular_grid(ar_shape, n_points):
     """Find `n_points` regularly spaced along `ar_shape`.
     
@@ -36,7 +37,7 @@ def regular_grid(ar_shape, n_points):
                                             (1.0 / (ndim - dim - 1)))
             if (sorted_dims >= stepsizes).all():
                 break
-    starts = np.floor(stepsizes/2)
+    starts = stepsizes // 2
     stepsizes = np.round(stepsizes)
     slices = [slice(start, None, step) for 
                                         start, step in zip(starts, stepsizes)]

From 441f4029f67db8733789f1d2e7f12e3c306b5055 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 14 May 2013 17:57:15 +1000
Subject: [PATCH 14/39] PEP8

---
 skimage/util/regular_grid.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/util/regular_grid.py b/skimage/util/regular_grid.py
index ab1b12dc..d1cb5428 100644
--- a/skimage/util/regular_grid.py
+++ b/skimage/util/regular_grid.py
@@ -34,12 +34,12 @@ def regular_grid(ar_shape, n_points):
             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)))
+                                 (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)]
+              start, step in zip(starts, stepsizes)]
     slices = [slices[i] for i in unsort_dim_idxs]
     return slices

From 99f905bc547319d58a34c97911c7e27cadf9ac29 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Wed, 15 May 2013 17:18:35 +1000
Subject: [PATCH 15/39] Always convert image to float in [0, 1]

---
 skimage/segmentation/slic.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/slic.py b/skimage/segmentation/slic.py
index d039d5ed..49f51773 100644
--- a/skimage/segmentation/slic.py
+++ b/skimage/segmentation/slic.py
@@ -50,8 +50,10 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
 
     Notes
     -----
-    The image is optionally smoothed using a Gaussian kernel prior to
-    segmentation.
+    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.
 
     References
     ----------
@@ -72,6 +74,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
             (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 and is_rgb(image):
@@ -80,7 +83,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     if not isinstance(sigma, coll.Iterable):
         sigma = np.array([sigma, sigma, sigma, 0])
     if (sigma > 0).any():
-        image = ndimage.gaussian_filter(img_as_float(image), sigma)
+        image = ndimage.gaussian_filter(image, sigma)
     if convert2lab:
         image = rgb2lab(image)
 

From e757b9ae06fff23615971c10e03ed52f2e87c83d Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 20 May 2013 14:46:11 +1000
Subject: [PATCH 16/39] Rename slic.py to avoid name conflicts

---
 skimage/segmentation/__init__.py                      | 2 +-
 skimage/segmentation/{slic.py => slic_superpixels.py} | 0
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename skimage/segmentation/{slic.py => slic_superpixels.py} (100%)

diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py
index ae1bf074..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.py b/skimage/segmentation/slic_superpixels.py
similarity index 100%
rename from skimage/segmentation/slic.py
rename to skimage/segmentation/slic_superpixels.py

From e9c5f666aed6566c7eba2ff04291590161fe6995 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 20 May 2013 16:30:55 +1000
Subject: [PATCH 17/39] Use educated guesses for presence of color channels in
 SLIC

---
 skimage/segmentation/slic_superpixels.py | 28 +++++++++++++++++++-----
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 49f51773..43639f16 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -1,14 +1,15 @@
 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, is_rgb
+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=True, convert2lab=True):
+         multichannel=None, convert2lab=True):
     """Segments image using k-means clustering in Color-(x,y) space.
 
     Parameters
@@ -26,9 +27,11 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     sigma : float, optional (default: 1)
         Width of Gaussian smoothing kernel for preprocessing. Zero means no
         smoothing.
-    multichannel : bool, optional (default: True)
+    multichannel : bool, optional (default: None)
         Whether the last axis of the image is to be interpreted as multiple
-        channels. Only 3 channels are supported.
+        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
@@ -46,7 +49,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
             - 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.
+            the image is not 3, OR
 
     Notes
     -----
@@ -55,6 +58,10 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
 
     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,
@@ -70,6 +77,15 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     >>> # 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 == image.ndim + 1)
     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)):
@@ -77,7 +93,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     image = img_as_float(image)
     if not multichannel:
         image = gray2rgb(image)
-    if image.ndim == 3 and is_rgb(image):
+    elif image.ndim == 3:
         # See 2D RGB image as 3D RGB image with Z = 1
         image = image[np.newaxis, ...]
     if not isinstance(sigma, coll.Iterable):

From 1c116a9905beddc188c0bef2a25fc0bdbe42872b Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 20 May 2013 16:42:06 +1000
Subject: [PATCH 18/39] Add function to guesstimate number of spatial
 dimensions

---
 skimage/color/__init__.py  |  2 ++
 skimage/color/colorconv.py | 31 +++++++++++++++++++++++++++++++
 2 files changed, 33 insertions(+)

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 540c5a75..5225f7ea 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.

From f0d586a32f5b902c66b9c8907a5039fde9538597 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 20 May 2013 16:43:20 +1000
Subject: [PATCH 19/39] Revert change to is_rgb

---
 skimage/color/colorconv.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index 5225f7ea..a73f5f19 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -90,7 +90,7 @@ def is_rgb(image):
         Input image.
 
     """
-    return (image.ndim in (3, 4) and image.shape[-1] in (3, 4))
+    return (image.ndim == 3 and image.shape[2] in (3, 4))
 
 
 @deprecated()

From 5a43e0306100fc8f5109000f6708a265853bf972 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 20 May 2013 16:55:07 +1000
Subject: [PATCH 20/39] Bug fix: add z dimension regardless of multichannel

---
 skimage/segmentation/slic_superpixels.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 43639f16..35e00a34 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -93,7 +93,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     image = img_as_float(image)
     if not multichannel:
         image = gray2rgb(image)
-    elif image.ndim == 3:
+    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):

From 55fe6ce2d685e25cfc9e0894c517bccd6b7bf4fe Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 20 May 2013 17:05:59 +1000
Subject: [PATCH 21/39] Bug fix: spatial vs image dimension comparison reversed

---
 skimage/segmentation/slic_superpixels.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 35e00a34..1c9b8510 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -85,7 +85,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
         warnings.warn(RuntimeWarning(msg))
         multichannel = True
     elif multichannel is None:
-        multichannel = (spatial_dims == image.ndim + 1)
+        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)):

From 69fb3fb7ba12055a0378330b26a3186e6638c615 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Thu, 30 May 2013 13:37:59 +1000
Subject: [PATCH 22/39] Make means array contiguous for numpy <= 1.6.1

---
 skimage/segmentation/slic_superpixels.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 1c9b8510..a4e949fd 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -117,6 +117,7 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
     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

From 8c2011d4f92858279d995373a0d6d2a6aa26312e Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 4 Jun 2013 14:25:40 +1000
Subject: [PATCH 23/39] Add UTF8 coding declaration to slic_superpixels.py

---
 skimage/segmentation/slic_superpixels.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index a4e949fd..2caf8ade 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -1,3 +1,5 @@
+# coding=utf-8
+
 import collections as coll
 import numpy as np
 from scipy import ndimage

From cba821d5e3128e0abb8f8fadbeb395531d6390a3 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Wed, 12 Jun 2013 23:14:22 +1000
Subject: [PATCH 24/39] Suppress warning of ambiguous array dim in test

---
 skimage/segmentation/tests/test_slic.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 1018969a..d0539cd2 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -1,4 +1,5 @@
 import itertools as it
+import warnings
 import numpy as np
 from numpy.testing import assert_equal, assert_array_equal
 from skimage.segmentation import slic
@@ -13,7 +14,9 @@ def test_color_2d():
     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)

From b1b70631bdecbfa7240a22aaa202ad88953baf78 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 24 Jun 2013 22:44:50 -0400
Subject: [PATCH 25/39] Initial attempt at updating SLIC for memoryviews

---
 skimage/segmentation/_slic.pyx | 44 +++++++++++-----------------------
 1 file changed, 14 insertions(+), 30 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index d668f084..99d8d51b 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -13,10 +13,10 @@ from ..util import img_as_float, regular_grid
 from ..color import rgb2lab, gray2rgb
 
 
-def _slic_cython(cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx,
-                 cnp.ndarray[dtype=cnp.intp_t, ndim=3] nearest_mean,
-                 cnp.ndarray[dtype=cnp.float_t, ndim=3] distance,
-                 cnp.ndarray[dtype=cnp.float_t, ndim=2] means,
+def _slic_cython(double[:, :, :, ::1] image_zyx,
+                 int[:, :, ::1] nearest_mean,
+                 double[:, :, ::1] distance,
+                 double[:, ::1] means,
                  float ratio, int max_iter, int n_segments):
     """Helper function for SLIC segmentation."""
 
@@ -30,54 +30,38 @@ def _slic_cython(cnp.ndarray[dtype=cnp.float_t, ndim=4] image_zyx,
     slices = regular_grid((depth, height, width), n_segments)
     step_z, step_y, step_x = [int(s.step) for s in slices]
 
-    cdef cnp.float_t* current_mean
-    cdef cnp.float_t* mean_entry
     n_means = means.shape[0]
     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.float_t* image_p = <cnp.float_t*> image_zyx.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
         # assign pixels to means
         for k in range(n_means):
             # compute windows:
-            z_min = int(max(current_mean[0] - 2 * step_z, 0))
-            z_max = int(min(current_mean[0] + 2 * step_z, depth))
-            y_min = int(max(current_mean[1] - 2 * step_y, 0))
-            y_max = int(min(current_mean[1] + 2 * step_y, height))
-            x_min = int(max(current_mean[2] - 2 * step_x, 0))
-            x_max = int(min(current_mean[2] + 2 * step_x, width))
+            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):
-                    current_pixel = \
-                            &image_p[6 * ((z * height + y) * width + x_min)]
-                    current_distance = \
-                            &distance_p[(z * height + y) * width + x_min]
                     for x in range(x_min, x_max):
-                        mean_entry = current_mean
                         dist_mean = 0
                         for c in range(6):
                             # you would think the compiler can optimize the
                             # squaring itself. mine can't (with O2)
-                            tmp = current_pixel[0] - mean_entry[0]
+                            tmp = image_zyx[z, y, x, c] - means[k, c]
                             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:
+                        if distance[z, y, x] - dist_mean > 1e-10:
                             nearest_mean[z, y, x] = k
-                            current_distance[0] = dist_mean
-                            changes += 1
-                        current_distance += 1
-            current_mean += 6
+                            distance[z, y, x] = dist_mean
+                            changes = 1
         if changes == 0:
             break
         # recompute means:

From 6e64515ea9060cea77a12cde7a0923ad56d4a206 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 24 Jun 2013 23:41:53 -0400
Subject: [PATCH 26/39] Fix build errors with memoryviews

---
 skimage/segmentation/_slic.pyx | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 99d8d51b..7f0fa8fb 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -14,7 +14,7 @@ from ..color import rgb2lab, gray2rgb
 
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
-                 int[:, :, ::1] nearest_mean,
+                 long[:, :, ::1] nearest_mean,
                  double[:, :, ::1] distance,
                  double[:, ::1] means,
                  float ratio, int max_iter, int n_segments):
@@ -22,8 +22,8 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
 
     # initialize on grid:
     cdef Py_ssize_t depth, height, width
-    shape = image_zyx.shape
-    depth, height, width = shape[0], shape[1], shape[2]
+    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_z, step_y, step_x
     grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
@@ -36,8 +36,10 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     cdef double dist_mean
 
     cdef double tmp
+    #cdef long[::1] nearest_mean_ravel
+    #cdef double[::1] image_zyx_ravel_j
     for i in range(max_iter):
-        distance.fill(np.inf)
+        distance[:, :, :] = np.inf
         changes = 0
         # assign pixels to means
         for k in range(n_means):
@@ -65,9 +67,13 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         if changes == 0:
             break
         # recompute means:
-        means_list = [np.bincount(nearest_mean.ravel(),
-                      image_zyx[:, :, :, j].ravel()) for j in range(6)]
-        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.asarray(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

From 524255c0a272c69d15d7703a498e9df147a2244b Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 1 Jul 2013 18:34:23 +0200
Subject: [PATCH 27/39] Improve _slic.pyx doc, bug fixes, debug print

---
 skimage/segmentation/_slic.pyx | 31 ++++++++++++++++++++++++++++---
 1 file changed, 28 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 7f0fa8fb..6621d4bf 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -18,7 +18,32 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                  double[:, :, ::1] distance,
                  double[:, ::1] means,
                  float ratio, int max_iter, int n_segments):
-    """Helper function for SLIC segmentation."""
+    """Helper function for SLIC segmentation.
+    
+    Parameters
+    ----------
+    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
+    -------
+    nearest_mean : 3D np.ndarray of long, shape (Z, Y, X)
+        The label field/superpixels found by SLIC.
+    """
 
     # initialize on grid:
     cdef Py_ssize_t depth, height, width
@@ -39,7 +64,6 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     #cdef long[::1] nearest_mean_ravel
     #cdef double[::1] image_zyx_ravel_j
     for i in range(max_iter):
-        distance[:, :, :] = np.inf
         changes = 0
         # assign pixels to means
         for k in range(n_means):
@@ -70,10 +94,11 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         nearest_mean_ravel = np.asarray(nearest_mean).ravel()
         means_list = []
         for j in range(6):
-            image_zyx_ravel = np.asarray(image_zyx[:, :, :, j]).ravel()
+            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")
+    print np.asarray(nearest_mean)
     return nearest_mean

From a33861896784b1c19d9d89050639f913eba65f49 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 1 Jul 2013 21:31:06 +0200
Subject: [PATCH 28/39] Remove diagnostic print

---
 skimage/segmentation/_slic.pyx | 1 -
 1 file changed, 1 deletion(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 6621d4bf..c7517dc3 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -100,5 +100,4 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         in_mean = np.bincount(nearest_mean_ravel)
         in_mean[in_mean == 0] = 1
         means = (np.vstack(means_list) / in_mean).T.copy("C")
-    print np.asarray(nearest_mean)
     return nearest_mean

From 578876eff89e9e5eb35135e9aead648a8624244b Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 1 Jul 2013 22:06:03 +0200
Subject: [PATCH 29/39] Bug fix: correctly initialize distance in slic

---
 skimage/segmentation/_slic.pyx | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index c7517dc3..337be882 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -65,6 +65,10 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     #cdef double[::1] image_zyx_ravel_j
     for i in range(max_iter):
         changes = 0
+        for z in range(depth):
+            for y in range(height):
+                for x in range(width):
+                    distance[z, y, x] = np.inf
         # assign pixels to means
         for k in range(n_means):
             # compute windows:
@@ -100,4 +104,4 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         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)

From cd197527fb39796520dd22e1f8df97b17a4810b3 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 2 Jul 2013 13:16:17 +0200
Subject: [PATCH 30/39] Speed up initialising 'distance' in SLIC

---
 skimage/segmentation/_slic.pyx | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 337be882..654c5cff 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -51,7 +51,6 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                             image_zyx.shape[2])
     # approximate grid size for desired n_segments
     cdef Py_ssize_t step_z, step_y, step_x
-    grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
     slices = regular_grid((depth, height, width), n_segments)
     step_z, step_y, step_x = [int(s.step) for s in slices]
 
@@ -65,10 +64,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     #cdef double[::1] image_zyx_ravel_j
     for i in range(max_iter):
         changes = 0
-        for z in range(depth):
-            for y in range(height):
-                for x in range(width):
-                    distance[z, y, x] = np.inf
+        distance[:, :, :] = np.inf
         # assign pixels to means
         for k in range(n_means):
             # compute windows:

From 7f1b9ea4c6f9f5ac1c0bff2d4d73aa0dd2387759 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 2 Jul 2013 19:42:24 +0200
Subject: [PATCH 31/39] Fix automerge bug

---
 skimage/util/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index 474b2966..b196af34 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -23,5 +23,5 @@ __all__ = ['img_as_float',
            'view_as_blocks',
            'view_as_windows',
            'pad',
-           'random_noise']
+           'random_noise',
            'regular_grid']

From dab831b72b59498de1cd2e30755f41d009826f50 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Tue, 2 Jul 2013 20:13:09 +0200
Subject: [PATCH 32/39] Bug fix: map returns iterator in Py3k

---
 skimage/color/colorconv.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index a73f5f19..d2b20316 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -739,7 +739,7 @@ def xyz2lab(xyz):
     a = 500.0 * (x - y)
     b = 200.0 * (y - z)
 
-    return np.concatenate(map(lambda x: x[..., np.newaxis], [L, a, b]), -1)
+    return np.concatenate([x[..., np.newaxis] for x in [L, a, b]], axis=-1)
 
 
 def lab2xyz(lab):

From df2233cc6ae9ee8dac800a204c656e5fa2efddc2 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Wed, 3 Jul 2013 00:49:36 +0200
Subject: [PATCH 33/39] Remove old commented-out code

---
 skimage/segmentation/_slic.pyx | 2 --
 1 file changed, 2 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 654c5cff..f6c6788c 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -60,8 +60,6 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     cdef double dist_mean
 
     cdef double tmp
-    #cdef long[::1] nearest_mean_ravel
-    #cdef double[::1] image_zyx_ravel_j
     for i in range(max_iter):
         changes = 0
         distance[:, :, :] = np.inf

From 48d78c09d94904047c87744a639b975646c11f24 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 8 Jul 2013 14:41:42 +0200
Subject: [PATCH 34/39] Improve doc for regular_grid.py

---
 skimage/util/regular_grid.py | 29 ++++++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/skimage/util/regular_grid.py b/skimage/util/regular_grid.py
index d1cb5428..e304be20 100644
--- a/skimage/util/regular_grid.py
+++ b/skimage/util/regular_grid.py
@@ -5,7 +5,10 @@ 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.
+    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
     ----------
@@ -20,6 +23,30 @@ def regular_grid(ar_shape, n_points):
     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)

From e80dd18d7e7dbdfeaad63b384724d5b0fa816e4e Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Fri, 12 Jul 2013 18:26:11 +0200
Subject: [PATCH 35/39] Add tests for 'regular_grid()'

---
 skimage/util/tests/test_regular_grid.py | 33 +++++++++++++++++++++++++
 1 file changed, 33 insertions(+)
 create mode 100644 skimage/util/tests/test_regular_grid.py

diff --git a/skimage/util/tests/test_regular_grid.py b/skimage/util/tests/test_regular_grid.py
new file mode 100644
index 00000000..1216b29d
--- /dev/null
+++ b/skimage/util/tests/test_regular_grid.py
@@ -0,0 +1,33 @@
+import numpy as np
+from nose.tools import raises
+from numpy.testing import assert_equal
+from skimage.util.regular_grid import regular_grid
+
+
+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()

From f0ecc4d00de22f7e05421cdd3eb3f525e7813eb6 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Fri, 12 Jul 2013 18:28:11 +0200
Subject: [PATCH 36/39] Rename regular_grid.py

---
 skimage/util/__init__.py                           | 2 +-
 skimage/util/{regular_grid.py => _regular_grid.py} | 0
 skimage/util/tests/test_regular_grid.py            | 2 +-
 3 files changed, 2 insertions(+), 2 deletions(-)
 rename skimage/util/{regular_grid.py => _regular_grid.py} (100%)

diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index b196af34..7afcf54a 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -11,7 +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
+from ._regular_grid import regular_grid
 
 
 __all__ = ['img_as_float',
diff --git a/skimage/util/regular_grid.py b/skimage/util/_regular_grid.py
similarity index 100%
rename from skimage/util/regular_grid.py
rename to skimage/util/_regular_grid.py
diff --git a/skimage/util/tests/test_regular_grid.py b/skimage/util/tests/test_regular_grid.py
index 1216b29d..2808d99c 100644
--- a/skimage/util/tests/test_regular_grid.py
+++ b/skimage/util/tests/test_regular_grid.py
@@ -1,7 +1,7 @@
 import numpy as np
 from nose.tools import raises
 from numpy.testing import assert_equal
-from skimage.util.regular_grid import regular_grid
+from skimage.util import regular_grid
 
 
 def test_regular_grid_2d_8():

From 6221ce6e3b86c8bb63837673df412d8d542a0f31 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Sat, 13 Jul 2013 01:27:17 +0200
Subject: [PATCH 37/39] Add test coverage for guess_spatial_dimensions()

---
 skimage/color/tests/test_colorconv.py | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/skimage/color/tests/test_colorconv.py b/skimage/color/tests/test_colorconv.py
index 05ab6915..6b06718b 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(guess_spatial_dimensions(im5), ValueError)
+
+
 class TestColorconv(TestCase):
 
     img_rgb = imread(os.path.join(data_dir, 'color.png'))

From 858f7411d722dc2ccf08f794e2b3bb0b2dcd35b0 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Sat, 13 Jul 2013 01:34:51 +0200
Subject: [PATCH 38/39] Complete test_regular_grid coverage

---
 skimage/util/tests/test_regular_grid.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/skimage/util/tests/test_regular_grid.py b/skimage/util/tests/test_regular_grid.py
index 2808d99c..61736a76 100644
--- a/skimage/util/tests/test_regular_grid.py
+++ b/skimage/util/tests/test_regular_grid.py
@@ -1,9 +1,16 @@
 import numpy as np
-from nose.tools import raises
 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)

From 8736a6650462dbeca176c19364917ce89a870ae2 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Sat, 13 Jul 2013 09:49:38 +1000
Subject: [PATCH 39/39] Fix typo using assert_raises in test_colorconv.py

---
 skimage/color/tests/test_colorconv.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/color/tests/test_colorconv.py b/skimage/color/tests/test_colorconv.py
index 6b06718b..4fdaa4c8 100644
--- a/skimage/color/tests/test_colorconv.py
+++ b/skimage/color/tests/test_colorconv.py
@@ -52,7 +52,7 @@ def test_guess_spatial_dimensions():
     assert_equal(guess_spatial_dimensions(im2), 3)
     assert_equal(guess_spatial_dimensions(im3), None)
     assert_equal(guess_spatial_dimensions(im4), 3)
-    assert_raises(guess_spatial_dimensions(im5), ValueError)
+    assert_raises(ValueError, guess_spatial_dimensions, im5)
 
 
 class TestColorconv(TestCase):