From fd729a4e30f87fc6b35ad5e81f241c7065f87683 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 17:28:48 +0200
Subject: [PATCH 01/29] Improve SLIC

---
 skimage/segmentation/_slic.pyx           | 80 ++++++++++++++----------
 skimage/segmentation/slic_superpixels.py | 41 +++++++-----
 skimage/segmentation/tests/test_slic.py  |  2 +-
 3 files changed, 73 insertions(+), 50 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 5df747ad..13e05bc2 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -2,31 +2,32 @@
 #cython: boundscheck=False
 #cython: nonecheck=False
 #cython: wraparound=False
-import numpy as np
-from scipy import ndimage
+from libc.float cimport DBL_MAX
 
-from skimage.util import img_as_float, regular_grid
-from skimage.color import rgb2lab, gray2rgb
+import numpy as np
+cimport numpy as cnp
+
+from skimage.util import regular_grid
 
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
                  Py_ssize_t[:, :, ::1] nearest_mean,
                  double[:, :, ::1] distance,
-                 double[:, ::1] means,
+                 double[:, ::1] clusters,
                  Py_ssize_t max_iter, Py_ssize_t n_segments):
     """Helper function for SLIC segmentation.
 
     Parameters
     ----------
-    image_zyx : 4D array of double, shape (Z, Y, X, 6)
+    image_zyx : 4D array of double, shape (Z, Y, X, C)
         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
+        `array([i, j, k, c])`, depending
         on the colorspace.
     nearest_mean : 3D array of int, shape (Z, Y, X)
         The (initially empty) label field.
     distance : 3D array of double, shape (Z, Y, X)
         The (initially infinity) array of distances to the nearest centroid.
-    means : 2D array of double, shape (n_segments, 6)
+    clusters : 2D array of double, shape (n_segments, 6)
         The centroids obtained by SLIC.
     max_iter : int
         The maximum number of k-means iterations.
@@ -43,36 +44,43 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     cdef Py_ssize_t depth, height, width
     depth, height, width = (image_zyx.shape[0], image_zyx.shape[1],
                             image_zyx.shape[2])
+
+    cdef Py_ssize_t n_features = clusters.shape[1]
+    cdef Py_ssize_t n_clusters = clusters.shape[0]
+
     # approximate grid size for desired n_segments
     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]
 
-    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 double tmp
+
+    cdef Py_ssize_t[:] n_cluster_elems = np.zeros(n_clusters, dtype=np.intp)
+
     for i in range(max_iter):
         changes = 0
-        distance[:, :, :] = np.inf
-        # assign pixels to means
-        for k in range(n_means):
+        distance[:, :, :] = DBL_MAX
+
+        # assign pixels to clusters
+        for k in range(n_clusters):
             # compute windows:
-            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))
+            z_min = int(max(clusters[k, 0] - 2 * step_z, 0))
+            z_max = int(min(clusters[k, 0] + 2 * step_z, depth))
+            y_min = int(max(clusters[k, 1] - 2 * step_y, 0))
+            y_max = int(min(clusters[k, 1] + 2 * step_y, height))
+            x_min = int(max(clusters[k, 2] - 2 * step_x, 0))
+            x_max = int(min(clusters[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):
+                        for c in range(n_features):
                             # 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]
+                            tmp = image_zyx[z, y, x, c] - clusters[k, c]
                             dist_mean += tmp * tmp
                         if distance[z, y, x] > dist_mean:
                             nearest_mean[z, y, x] = k
@@ -80,15 +88,23 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                             changes = 1
         if changes == 0:
             break
-        # recompute means:
-        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 np.ascontiguousarray(nearest_mean)
+
+        # recompute clusters
+
+        # sum features for all clusters
+        n_cluster_elems[:] = 0
+        clusters[:, :] = 0
+        for z in range(depth):
+            for y in range(height):
+                for x in range(width):
+                    k = nearest_mean[z, y, x]
+                    n_cluster_elems[k] += 1
+                    for c in range(n_features):
+                        clusters[k, c] += image_zyx[z, y, x, c]
+
+        # divide by number of elements per cluster to obtain mean
+        for k in range(n_clusters):
+            for c in range(n_features):
+                clusters[k, c] /= n_cluster_elems[k]
+
+    return np.asarray(nearest_mean)
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index b905f199..28f6ef8f 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -21,7 +21,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
         (see `multichannel` parameter).
     n_segments : int, optional (default: 100)
         The (approximate) number of labels in the segmented output image.
-    compactness: float, optional (default: 10)
+    compactness : float, optional (default: 10)
         Balances color-space proximity and image-space proximity. Higher
         values give more weight to image-space. As `compactness` tends to
         infinity, superpixel shapes become square/cubic.
@@ -37,14 +37,14 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
         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
+        segmentation. For this purpose, the input is assumed to be RGB. Highly
         recommended.
     ratio : float, optional
         Synonym for `compactness`. This keyword is deprecated.
 
     Returns
     -------
-    segment_mask : (width, height) ndarray
+    segment_mask : (width, height, depth) array
         Integer mask indicating segment labels.
 
     Raises
@@ -99,34 +99,41 @@ def slic(image, n_segments=100, compactness=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)
+    image = np.atleast_3d(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:
+
+    if image.shape[3] == 3 and convert2lab:
         image = rgb2lab(image)
 
-    # initialize on grid:
+    # 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]
+    clusters_z = grid_z[slices]
+    clusters_y = grid_y[slices]
+    clusters_x = grid_x[slices]
+
+    clusters_color = np.zeros(clusters_z.shape + (image.shape[3],))
+    clusters = np.concatenate([clusters_z[..., np.newaxis],
+                               clusters_y[..., np.newaxis],
+                               clusters_x[..., np.newaxis],
+                               clusters_color
+                              ], axis=-1).reshape(-1, 3 + image.shape[3])
+    clusters = np.ascontiguousarray(clusters)
 
-    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 = float(max((step_z, step_y, step_x))) / compactness
@@ -134,9 +141,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
                                 grid_y[..., np.newaxis],
                                 grid_x[..., np.newaxis],
                                 image * ratio], axis=-1).copy("C")
-    nearest_mean = np.zeros((depth, height, width), dtype=np.intp)
+    nearest_cluster = np.empty((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,
+    segment_map = _slic_cython(image_zyx, nearest_cluster, distance, clusters,
                                max_iter, n_segments)
     if segment_map.shape[0] == 1:
         segment_map = segment_map[0]
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 9e6a39e2..6ecd928d 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -28,7 +28,7 @@ def test_color_2d():
 
 def test_gray_2d():
     rnd = np.random.RandomState(0)
-    img = np.zeros((20, 21))
+    img = np.zeros((20, 20))
     img[:10, :10] = 0.33
     img[10:, :10] = 0.67
     img[10:, 10:] = 1.00

From 6b5556b27969fd051b8c9046fb8b8c29f8dd813e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 18:25:18 +0200
Subject: [PATCH 02/29] Reduce memory footprint

---
 skimage/segmentation/_slic.pyx           | 67 +++++++++++++++---------
 skimage/segmentation/slic_superpixels.py | 11 ++--
 2 files changed, 48 insertions(+), 30 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 13e05bc2..86654b11 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -11,7 +11,7 @@ from skimage.util import regular_grid
 
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
-                 Py_ssize_t[:, :, ::1] nearest_mean,
+                 Py_ssize_t[:, :, ::1] nearest_clusters,
                  double[:, :, ::1] distance,
                  double[:, ::1] clusters,
                  Py_ssize_t max_iter, Py_ssize_t n_segments):
@@ -21,13 +21,12 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     ----------
     image_zyx : 4D array of double, shape (Z, Y, X, C)
         The image with embedded coordinates, that is, `image_zyx[i, j, k]` is
-        `array([i, j, k, c])`, depending
-        on the colorspace.
-    nearest_mean : 3D array of int, shape (Z, Y, X)
+        `array([i, j, k, c])`, depending on the colorspace.
+    nearest_clusters : 3D array of int, shape (Z, Y, X)
         The (initially empty) label field.
     distance : 3D array of double, shape (Z, Y, X)
         The (initially infinity) array of distances to the nearest centroid.
-    clusters : 2D array of double, shape (n_segments, 6)
+    clusters : 2D array of double, shape (n_segments, 3 + C)
         The centroids obtained by SLIC.
     max_iter : int
         The maximum number of k-means iterations.
@@ -36,7 +35,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
 
     Returns
     -------
-    nearest_mean : 3D array of int, shape (Z, Y, X)
+    nearest_clusters : 3D array of int, shape (Z, Y, X)
         The label field/superpixels found by SLIC.
     """
 
@@ -45,48 +44,63 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     depth, height, width = (image_zyx.shape[0], image_zyx.shape[1],
                             image_zyx.shape[2])
 
-    cdef Py_ssize_t n_features = clusters.shape[1]
     cdef Py_ssize_t n_clusters = clusters.shape[0]
+    # number of features [X, Y, Z, ...]
+    cdef Py_ssize_t n_features = clusters.shape[1]
 
     # approximate grid size for desired n_segments
     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]
 
-    cdef Py_ssize_t i, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max, \
-            changes
+    cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, \
+                    z_max
     cdef double dist_mean
-    cdef double tmp
+
+    cdef char change
+
+    cdef double cx, cy, cz, dy, dz
 
     cdef Py_ssize_t[:] n_cluster_elems = np.zeros(n_clusters, dtype=np.intp)
 
     for i in range(max_iter):
-        changes = 0
+        change = 0
         distance[:, :, :] = DBL_MAX
 
         # assign pixels to clusters
         for k in range(n_clusters):
-            # compute windows:
+
+            # compute windows
             z_min = int(max(clusters[k, 0] - 2 * step_z, 0))
             z_max = int(min(clusters[k, 0] + 2 * step_z, depth))
             y_min = int(max(clusters[k, 1] - 2 * step_y, 0))
             y_max = int(min(clusters[k, 1] + 2 * step_y, height))
             x_min = int(max(clusters[k, 2] - 2 * step_x, 0))
             x_max = int(min(clusters[k, 2] + 2 * step_x, width))
+
+            # cluster coordinate centers
+            cz = clusters[k, 0]
+            cy = clusters[k, 1]
+            cx = clusters[k, 2]
+
             for z in range(z_min, z_max):
+                dz = (cz - z) ** 2
                 for y in range(y_min, y_max):
+                    dy = (cy - y) ** 2
                     for x in range(x_min, x_max):
-                        dist_mean = 0
-                        for c in range(n_features):
-                            # you would think the compiler can optimize the
-                            # squaring itself. mine can't (with O2)
-                            tmp = image_zyx[z, y, x, c] - clusters[k, c]
-                            dist_mean += tmp * tmp
+                        # you would think the compiler can optimize the
+                        # squaring itself. mine can't (with O2)
+                        dist_mean = dz + dy + (cx - x) ** 2
+                        for c in range(3, n_features):
+                            dist_mean += (image_zyx[z, y, x, c - 3]
+                                          - clusters[k, c]) ** 2
                         if distance[z, y, x] > dist_mean:
-                            nearest_mean[z, y, x] = k
+                            nearest_clusters[z, y, x] = k
                             distance[z, y, x] = dist_mean
-                            changes = 1
-        if changes == 0:
+                            change = 1
+
+        # stop if no pixel changed its cluster
+        if change == 0:
             break
 
         # recompute clusters
@@ -97,14 +111,17 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         for z in range(depth):
             for y in range(height):
                 for x in range(width):
-                    k = nearest_mean[z, y, x]
+                    k = nearest_clusters[z, y, x]
                     n_cluster_elems[k] += 1
-                    for c in range(n_features):
-                        clusters[k, c] += image_zyx[z, y, x, c]
+                    clusters[k, 0] += z
+                    clusters[k, 1] += y
+                    clusters[k, 2] += x
+                    for c in range(3, n_features):
+                        clusters[k, c] += image_zyx[z, y, x, c - 3]
 
         # divide by number of elements per cluster to obtain mean
         for k in range(n_clusters):
             for c in range(n_features):
                 clusters[k, c] /= n_cluster_elems[k]
 
-    return np.asarray(nearest_mean)
+    return np.asarray(nearest_clusters)
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 28f6ef8f..1db2754a 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -137,14 +137,15 @@ def slic(image, n_segments=100, compactness=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 = float(max((step_z, step_y, step_x))) / compactness
-    image_zyx = np.concatenate([grid_z[..., np.newaxis],
-                                grid_y[..., np.newaxis],
-                                grid_x[..., np.newaxis],
-                                image * ratio], axis=-1).copy("C")
+    image = image * ratio
+
     nearest_cluster = np.empty((depth, height, width), dtype=np.intp)
     distance = np.empty((depth, height, width), dtype=np.float)
-    segment_map = _slic_cython(image_zyx, nearest_cluster, distance, clusters,
+
+    segment_map = _slic_cython(image, nearest_cluster, distance, clusters,
                                max_iter, n_segments)
+
     if segment_map.shape[0] == 1:
         segment_map = segment_map[0]
+
     return segment_map

From 471b293a9e982dbea5b11e1ddeb1b5da6a265293 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 18:34:18 +0200
Subject: [PATCH 03/29] Remove unnecessary parameters and update doc string

---
 skimage/segmentation/_slic.pyx           | 26 ++++++++++--------------
 skimage/segmentation/slic_superpixels.py | 18 +++++++---------
 2 files changed, 18 insertions(+), 26 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 86654b11..6faaa69b 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -11,27 +11,18 @@ from skimage.util import regular_grid
 
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
-                 Py_ssize_t[:, :, ::1] nearest_clusters,
-                 double[:, :, ::1] distance,
                  double[:, ::1] clusters,
-                 Py_ssize_t max_iter, Py_ssize_t n_segments):
+                 Py_ssize_t max_iter):
     """Helper function for SLIC segmentation.
 
     Parameters
     ----------
     image_zyx : 4D array of double, shape (Z, Y, X, C)
-        The image with embedded coordinates, that is, `image_zyx[i, j, k]` is
-        `array([i, j, k, c])`, depending on the colorspace.
-    nearest_clusters : 3D array of int, shape (Z, Y, X)
-        The (initially empty) label field.
-    distance : 3D array of double, shape (Z, Y, X)
-        The (initially infinity) array of distances to the nearest centroid.
-    clusters : 2D array of double, shape (n_segments, 3 + C)
-        The centroids obtained by SLIC.
+        The input image.
+    clusters : 2D array of double, shape (N, 3 + C)
+        The initial centroids obtained by SLIC as [Z, Y, X, C...].
     max_iter : int
         The maximum number of k-means iterations.
-    n_segments : int
-        The approximate/desired number of segments.
 
     Returns
     -------
@@ -39,7 +30,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         The label field/superpixels found by SLIC.
     """
 
-    # initialize on grid:
+    # initialize on grid
     cdef Py_ssize_t depth, height, width
     depth, height, width = (image_zyx.shape[0], image_zyx.shape[1],
                             image_zyx.shape[2])
@@ -50,9 +41,14 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
 
     # approximate grid size for desired n_segments
     cdef Py_ssize_t step_z, step_y, step_x
-    slices = regular_grid((depth, height, width), n_segments)
+    slices = regular_grid((depth, height, width), n_clusters)
     step_z, step_y, step_x = [int(s.step) for s in slices]
 
+    cdef Py_ssize_t[:, :, ::1] nearest_clusters \
+        = np.empty((depth, height, width), dtype=np.intp)
+    cdef float[:, :, ::1] distance \
+        = np.empty((depth, height, width), dtype=np.float32)
+
     cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, \
                     z_max
     cdef double dist_mean
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 1db2754a..42d30593 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -44,7 +44,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
 
     Returns
     -------
-    segment_mask : (width, height, depth) array
+    labels : (width, height, depth) array
         Integer mask indicating segment labels.
 
     Raises
@@ -82,6 +82,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
     >>> # Increasing the ratio parameter yields more square regions
     >>> segments = slic(img, n_segments=100, ratio=20)
     """
+
     if ratio is not None:
         msg = 'Keyword `ratio` is deprecated. Use `compactness` instead.'
         warnings.warn(msg)
@@ -115,10 +116,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
     if image.shape[3] == 3 and convert2lab:
         image = rgb2lab(image)
 
-    # initialize on grid
     depth, height, width = image.shape[:3]
 
-    # approximate grid size for desired n_segments
+    # initialize cluster centroids for desired number of 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]
@@ -139,13 +139,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
     ratio = float(max((step_z, step_y, step_x))) / compactness
     image = image * ratio
 
-    nearest_cluster = np.empty((depth, height, width), dtype=np.intp)
-    distance = np.empty((depth, height, width), dtype=np.float)
+    labels = _slic_cython(image, clusters, max_iter)
 
-    segment_map = _slic_cython(image, nearest_cluster, distance, clusters,
-                               max_iter, n_segments)
+    if labels.shape[0] == 1:
+        labels = labels[0]
 
-    if segment_map.shape[0] == 1:
-        segment_map = segment_map[0]
-
-    return segment_map
+    return labels

From 84579a6c2c82894ef88e88e6048cb198243a8265 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 18:37:35 +0200
Subject: [PATCH 04/29] Remove legacy comment

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 6faaa69b..030bcb4e 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -84,8 +84,6 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                 for y in range(y_min, y_max):
                     dy = (cy - y) ** 2
                     for x in range(x_min, x_max):
-                        # you would think the compiler can optimize the
-                        # squaring itself. mine can't (with O2)
                         dist_mean = dz + dy + (cx - x) ** 2
                         for c in range(3, n_features):
                             dist_mean += (image_zyx[z, y, x, c - 3]

From bc7efb01e4fe7e0d6f9a07c205101adf01f92ed8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 19:38:11 +0200
Subject: [PATCH 05/29] Fix dtype bug

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 030bcb4e..4a57ee93 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -46,8 +46,8 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
 
     cdef Py_ssize_t[:, :, ::1] nearest_clusters \
         = np.empty((depth, height, width), dtype=np.intp)
-    cdef float[:, :, ::1] distance \
-        = np.empty((depth, height, width), dtype=np.float32)
+    cdef double[:, :, ::1] distance \
+        = np.empty((depth, height, width), dtype=np.double)
 
     cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, \
                     z_max

From bd388383a02327250dee00481c7b0eab66d193fd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 19:59:00 +0200
Subject: [PATCH 06/29] Fix bug in window extent determination

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 4a57ee93..c61455e5 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -66,19 +66,19 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         # assign pixels to clusters
         for k in range(n_clusters):
 
-            # compute windows
-            z_min = int(max(clusters[k, 0] - 2 * step_z, 0))
-            z_max = int(min(clusters[k, 0] + 2 * step_z, depth))
-            y_min = int(max(clusters[k, 1] - 2 * step_y, 0))
-            y_max = int(min(clusters[k, 1] + 2 * step_y, height))
-            x_min = int(max(clusters[k, 2] - 2 * step_x, 0))
-            x_max = int(min(clusters[k, 2] + 2 * step_x, width))
-
             # cluster coordinate centers
             cz = clusters[k, 0]
             cy = clusters[k, 1]
             cx = clusters[k, 2]
 
+            # compute windows
+            z_min = <Py_ssize_t>max(cz - 2 * step_z, 0)
+            z_max = <Py_ssize_t>min(cz + 2 * step_z + 1, depth)
+            y_min = <Py_ssize_t>max(cy - 2 * step_y, 0)
+            y_max = <Py_ssize_t>min(cy + 2 * step_y + 1, height)
+            x_min = <Py_ssize_t>max(cx - 2 * step_x, 0)
+            x_max = <Py_ssize_t>min(cx + 2 * step_x + 1, width)
+
             for z in range(z_min, z_max):
                 dz = (cz - z) ** 2
                 for y in range(y_min, y_max):

From f5cfbcfe97f17f5e9c6ce782e6869eaf96dc7016 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 31 Aug 2013 20:04:36 +0200
Subject: [PATCH 07/29] Reorder variable declarations

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index c61455e5..834c5181 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -32,8 +32,9 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
 
     # initialize on grid
     cdef Py_ssize_t depth, height, width
-    depth, height, width = (image_zyx.shape[0], image_zyx.shape[1],
-                            image_zyx.shape[2])
+    depth = image_zyx.shape[0]
+    height = image_zyx.shape[1]
+    width = image_zyx.shape[2]
 
     cdef Py_ssize_t n_clusters = clusters.shape[0]
     # number of features [X, Y, Z, ...]
@@ -48,17 +49,12 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         = np.empty((depth, height, width), dtype=np.intp)
     cdef double[:, :, ::1] distance \
         = np.empty((depth, height, width), dtype=np.double)
-
-    cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, \
-                    z_max
-    cdef double dist_mean
-
-    cdef char change
-
-    cdef double cx, cy, cz, dy, dz
-
     cdef Py_ssize_t[:] n_cluster_elems = np.zeros(n_clusters, dtype=np.intp)
 
+    cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max
+    cdef char change
+    cdef double dist_mean, cx, cy, cz, dy, dz
+
     for i in range(max_iter):
         change = 0
         distance[:, :, :] = DBL_MAX

From e5eea8e135264dacac28be5bf6d468f1f3cc2c57 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sun, 1 Sep 2013 16:15:00 +0200
Subject: [PATCH 08/29] Improve documentation of sigma

---
 skimage/segmentation/slic_superpixels.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 42d30593..8522b482 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -27,9 +27,10 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
         infinity, superpixel shapes become square/cubic.
     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.
+    sigma : float or array of floats, optional (default: 1)
+        Width of Gaussian smoothing kernel for pre-processing for each
+        dimension of the image. The same sigma is applied to each dimension in
+        case of a scalar value. 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

From eeddd9e35fabb22d399da7cc564f25b9276705a6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sun, 1 Sep 2013 16:44:59 +0200
Subject: [PATCH 09/29] Revert multichannel magic and improve parameter docs

---
 skimage/segmentation/slic_superpixels.py | 58 ++++++++++--------------
 skimage/segmentation/tests/test_slic.py  | 10 ++--
 2 files changed, 30 insertions(+), 38 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 8522b482..c239baf0 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -6,37 +6,34 @@ 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, compactness=10., max_iter=10, sigma=1,
-         multichannel=None, convert2lab=True, ratio=None):
-    """Segments image using k-means clustering in Color-(x,y) space.
+def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=1,
+         multichannel=True, convert2lab=True, ratio=None):
+    """Segments image using k-means clustering in Color-(x,y,z) space.
 
     Parameters
     ----------
-    image : (width, height [, depth] [, 3]) ndarray
-        Input image, which can be 2D or 3D, and grayscale or multi-channel
+    image : 2D, 3D or 4D ndarray
+        Input image, which can be 2D or 3D, and grayscale or multichannel
         (see `multichannel` parameter).
-    n_segments : int, optional (default: 100)
+    n_segments : int
         The (approximate) number of labels in the segmented output image.
-    compactness : float, optional (default: 10)
+    compactness : float
         Balances color-space proximity and image-space proximity. Higher
         values give more weight to image-space. As `compactness` tends to
         infinity, superpixel shapes become square/cubic.
-    max_iter : int, optional (default: 10)
+    max_iter : int
         Maximum number of iterations of k-means.
-    sigma : float or array of floats, optional (default: 1)
+    sigma : float or (3,) array of floats
         Width of Gaussian smoothing kernel for pre-processing for each
         dimension of the image. The same sigma is applied to each dimension in
         case of a scalar value. Zero means no smoothing.
-    multichannel : bool, optional (default: None)
+    multichannel : bool
         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)
+        channels or another spatial dimension.
+    convert2lab : bool
         Whether the input should be converted to Lab colorspace prior to
         segmentation. For this purpose, the input is assumed to be RGB. Highly
         recommended.
@@ -45,7 +42,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
 
     Returns
     -------
-    labels : (width, height, depth) array
+    labels : 2D or 3D array
         Integer mask indicating segment labels.
 
     Raises
@@ -88,33 +85,28 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1,
         msg = 'Keyword `ratio` is deprecated. Use `compactness` instead.'
         warnings.warn(msg)
         compactness = ratio
-    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)
     image = np.atleast_3d(image)
 
     if image.ndim == 3:
-        # See 2D RGB image as 3D RGB image with Z = 1
-        image = image[np.newaxis, ...]
+        if multichannel:
+            # Make 2D image 3D with depth = 1
+            image = image[np.newaxis, ...]
+        else:
+            # Add channel as single last dimension
+            image = image[..., np.newaxis]
 
     if not isinstance(sigma, coll.Iterable):
-        sigma = np.array([sigma, sigma, sigma, 0])
+        sigma = np.array([sigma, sigma, sigma])
     if (sigma > 0).any():
+        sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
 
-    if image.shape[3] == 3 and convert2lab:
+    if convert2lab:
+
+        if not multichannel or image.shape[3] != 3:
+            raise ValueError("Lab colorspace conversion requires a RGB image.")
         image = rgb2lab(image)
 
     depth, height, width = image.shape[:3]
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 6ecd928d..cb16e1b7 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -28,15 +28,15 @@ def test_color_2d():
 
 def test_gray_2d():
     rnd = np.random.RandomState(0)
-    img = np.zeros((20, 20))
+    img = np.zeros((20, 21))
     img[:10, :10] = 0.33
     img[10:, :10] = 0.67
     img[10:, 10:] = 1.00
     img += 0.0033 * rnd.normal(size=img.shape)
     img[img > 1] = 1
     img[img < 0] = 0
-    seg = slic(img, sigma=0, n_segments=4, compactness=20.0,
-               multichannel=False)
+    seg = slic(img, sigma=0, n_segments=4, compactness=1,
+               multichannel=False, convert2lab=False)
 
     assert_equal(len(np.unique(seg)), 4)
     assert_array_equal(seg[:10, :10], 0)
@@ -80,8 +80,8 @@ 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, compactness=20.0,
-               multichannel=False)
+    seg = slic(img, sigma=0, n_segments=8, compactness=1,
+               multichannel=False, convert2lab=False)
 
     assert_equal(len(np.unique(seg)), 8)
     for s, c in zip(slices, range(8)):

From a6d1b10e254912ab8a3dbb0c377842cbdd98d52a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sun, 1 Sep 2013 16:53:05 +0200
Subject: [PATCH 10/29] Use absolute imports

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

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index c239baf0..6e009120 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -5,8 +5,8 @@ import numpy as np
 from scipy import ndimage
 import warnings
 
-from ..util import img_as_float, regular_grid
-from ._slic import _slic_cython
+from skimage.util import img_as_float, regular_grid
+from skimage.segmentation._slic import _slic_cython
 
 
 def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=1,

From 80ed8751474d705706bb5469ad73108470c81b9a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sun, 1 Sep 2013 17:12:10 +0200
Subject: [PATCH 11/29] Add missing rgb2lab import

---
 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 6e009120..5a6d4506 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -7,6 +7,7 @@ import warnings
 
 from skimage.util import img_as_float, regular_grid
 from skimage.segmentation._slic import _slic_cython
+from skimage.color import rgb2lab
 
 
 def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=1,

From d9e64b43e6c3d3434ba49fef971a33b301436407 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Tue, 3 Sep 2013 08:40:18 +0200
Subject: [PATCH 12/29] Make image C-contiguous

---
 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 5a6d4506..4f10e75f 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -131,7 +131,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=20, 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 = float(max((step_z, step_y, step_x))) / compactness
-    image = image * ratio
+    image = np.ascontiguousarray(image * ratio)
 
     labels = _slic_cython(image, clusters, max_iter)
 

From 770e28d2bb6a37a7a32f9ed2a151a6a44b44e0a9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Tue, 3 Sep 2013 10:47:49 +0200
Subject: [PATCH 13/29] Rename clusters to segments

---
 skimage/segmentation/_slic.pyx           | 72 ++++++++++++------------
 skimage/segmentation/slic_superpixels.py | 20 +++----
 2 files changed, 46 insertions(+), 46 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 834c5181..ac45a5c4 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -11,7 +11,7 @@ from skimage.util import regular_grid
 
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
-                 double[:, ::1] clusters,
+                 double[:, ::1] segments,
                  Py_ssize_t max_iter):
     """Helper function for SLIC segmentation.
 
@@ -19,14 +19,14 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     ----------
     image_zyx : 4D array of double, shape (Z, Y, X, C)
         The input image.
-    clusters : 2D array of double, shape (N, 3 + C)
+    segments : 2D array of double, shape (N, 3 + C)
         The initial centroids obtained by SLIC as [Z, Y, X, C...].
     max_iter : int
         The maximum number of k-means iterations.
 
     Returns
     -------
-    nearest_clusters : 3D array of int, shape (Z, Y, X)
+    nearest_segments : 3D array of int, shape (Z, Y, X)
         The label field/superpixels found by SLIC.
     """
 
@@ -36,36 +36,36 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     height = image_zyx.shape[1]
     width = image_zyx.shape[2]
 
-    cdef Py_ssize_t n_clusters = clusters.shape[0]
+    cdef Py_ssize_t n_segments = segments.shape[0]
     # number of features [X, Y, Z, ...]
-    cdef Py_ssize_t n_features = clusters.shape[1]
+    cdef Py_ssize_t n_features = segments.shape[1]
 
     # approximate grid size for desired n_segments
     cdef Py_ssize_t step_z, step_y, step_x
-    slices = regular_grid((depth, height, width), n_clusters)
+    slices = regular_grid((depth, height, width), n_segments)
     step_z, step_y, step_x = [int(s.step) for s in slices]
 
-    cdef Py_ssize_t[:, :, ::1] nearest_clusters \
+    cdef Py_ssize_t[:, :, ::1] nearest_segments \
         = np.empty((depth, height, width), dtype=np.intp)
     cdef double[:, :, ::1] distance \
         = np.empty((depth, height, width), dtype=np.double)
-    cdef Py_ssize_t[:] n_cluster_elems = np.zeros(n_clusters, dtype=np.intp)
+    cdef Py_ssize_t[:] n_segment_elems = np.zeros(n_segments, dtype=np.intp)
 
     cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max
     cdef char change
-    cdef double dist_mean, cx, cy, cz, dy, dz
+    cdef double dist_center, cx, cy, cz, dy, dz
 
     for i in range(max_iter):
         change = 0
         distance[:, :, :] = DBL_MAX
 
-        # assign pixels to clusters
-        for k in range(n_clusters):
+        # assign pixels to segments
+        for k in range(n_segments):
 
-            # cluster coordinate centers
-            cz = clusters[k, 0]
-            cy = clusters[k, 1]
-            cx = clusters[k, 2]
+            # segment coordinate centers
+            cz = segments[k, 0]
+            cy = segments[k, 1]
+            cx = segments[k, 2]
 
             # compute windows
             z_min = <Py_ssize_t>max(cz - 2 * step_z, 0)
@@ -80,38 +80,38 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                 for y in range(y_min, y_max):
                     dy = (cy - y) ** 2
                     for x in range(x_min, x_max):
-                        dist_mean = dz + dy + (cx - x) ** 2
+                        dist_center = dz + dy + (cx - x) ** 2
                         for c in range(3, n_features):
-                            dist_mean += (image_zyx[z, y, x, c - 3]
-                                          - clusters[k, c]) ** 2
-                        if distance[z, y, x] > dist_mean:
-                            nearest_clusters[z, y, x] = k
-                            distance[z, y, x] = dist_mean
+                            dist_center += (image_zyx[z, y, x, c - 3]
+                                          - segments[k, c]) ** 2
+                        if distance[z, y, x] > dist_center:
+                            nearest_segments[z, y, x] = k
+                            distance[z, y, x] = dist_center
                             change = 1
 
-        # stop if no pixel changed its cluster
+        # stop if no pixel changed its segment
         if change == 0:
             break
 
-        # recompute clusters
+        # recompute segment centers
 
-        # sum features for all clusters
-        n_cluster_elems[:] = 0
-        clusters[:, :] = 0
+        # sum features for all segments
+        n_segment_elems[:] = 0
+        segments[:, :] = 0
         for z in range(depth):
             for y in range(height):
                 for x in range(width):
-                    k = nearest_clusters[z, y, x]
-                    n_cluster_elems[k] += 1
-                    clusters[k, 0] += z
-                    clusters[k, 1] += y
-                    clusters[k, 2] += x
+                    k = nearest_segments[z, y, x]
+                    n_segment_elems[k] += 1
+                    segments[k, 0] += z
+                    segments[k, 1] += y
+                    segments[k, 2] += x
                     for c in range(3, n_features):
-                        clusters[k, c] += image_zyx[z, y, x, c - 3]
+                        segments[k, c] += image_zyx[z, y, x, c - 3]
 
-        # divide by number of elements per cluster to obtain mean
-        for k in range(n_clusters):
+        # divide by number of elements per segment to obtain mean
+        for k in range(n_segments):
             for c in range(n_features):
-                clusters[k, c] /= n_cluster_elems[k]
+                segments[k, c] /= n_segment_elems[k]
 
-    return np.asarray(nearest_clusters)
+    return np.asarray(nearest_segments)
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 4f10e75f..c5470992 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -116,24 +116,24 @@ def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=1,
     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]
-    clusters_z = grid_z[slices]
-    clusters_y = grid_y[slices]
-    clusters_x = grid_x[slices]
+    segments_z = grid_z[slices]
+    segments_y = grid_y[slices]
+    segments_x = grid_x[slices]
 
-    clusters_color = np.zeros(clusters_z.shape + (image.shape[3],))
-    clusters = np.concatenate([clusters_z[..., np.newaxis],
-                               clusters_y[..., np.newaxis],
-                               clusters_x[..., np.newaxis],
-                               clusters_color
+    segments_color = np.zeros(segments_z.shape + (image.shape[3],))
+    segments = np.concatenate([segments_z[..., np.newaxis],
+                               segments_y[..., np.newaxis],
+                               segments_x[..., np.newaxis],
+                               segments_color
                               ], axis=-1).reshape(-1, 3 + image.shape[3])
-    clusters = np.ascontiguousarray(clusters)
+    segments = np.ascontiguousarray(segments)
 
     # we do the scaling of ratio in the same way as in the SLIC paper
     # so the values have the same meaning
     ratio = float(max((step_z, step_y, step_x))) / compactness
     image = np.ascontiguousarray(image * ratio)
 
-    labels = _slic_cython(image, clusters, max_iter)
+    labels = _slic_cython(image, segments, max_iter)
 
     if labels.shape[0] == 1:
         labels = labels[0]

From 8137c41e222c72941454a25fb0e214f6fc0c1ba9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Tue, 3 Sep 2013 10:48:53 +0200
Subject: [PATCH 14/29] Add optional description to parameters

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

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index c5470992..1adba964 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -19,22 +19,22 @@ def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=1,
     image : 2D, 3D or 4D ndarray
         Input image, which can be 2D or 3D, and grayscale or multichannel
         (see `multichannel` parameter).
-    n_segments : int
+    n_segments : int, optional
         The (approximate) number of labels in the segmented output image.
-    compactness : float
+    compactness : float, optional
         Balances color-space proximity and image-space proximity. Higher
         values give more weight to image-space. As `compactness` tends to
         infinity, superpixel shapes become square/cubic.
-    max_iter : int
+    max_iter : int, optional
         Maximum number of iterations of k-means.
-    sigma : float or (3,) array of floats
+    sigma : float or (3,) array of floats, optional
         Width of Gaussian smoothing kernel for pre-processing for each
         dimension of the image. The same sigma is applied to each dimension in
         case of a scalar value. Zero means no smoothing.
-    multichannel : bool
+    multichannel : bool, optional
         Whether the last axis of the image is to be interpreted as multiple
         channels or another spatial dimension.
-    convert2lab : bool
+    convert2lab : bool, optional
         Whether the input should be converted to Lab colorspace prior to
         segmentation. For this purpose, the input is assumed to be RGB. Highly
         recommended.

From 22ca0ae45729af1bb4a6bfe261f758b408aa46ce Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Tue, 3 Sep 2013 20:31:24 +0200
Subject: [PATCH 15/29] Fix indentation

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index ac45a5c4..d818a501 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -83,7 +83,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                         dist_center = dz + dy + (cx - x) ** 2
                         for c in range(3, n_features):
                             dist_center += (image_zyx[z, y, x, c - 3]
-                                          - segments[k, c]) ** 2
+                                            - segments[k, c]) ** 2
                         if distance[z, y, x] > dist_center:
                             nearest_segments[z, y, x] = k
                             distance[z, y, x] = dist_center

From 05aaeef5ac71d14f4d54eef2d9b8f5fb3c445346 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 14 Sep 2013 10:08:07 +0200
Subject: [PATCH 16/29] Change default value of sigma

---
 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 1adba964..4c4a7921 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -10,7 +10,7 @@ from skimage.segmentation._slic import _slic_cython
 from skimage.color import rgb2lab
 
 
-def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=1,
+def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=0,
          multichannel=True, convert2lab=True, ratio=None):
     """Segments image using k-means clustering in Color-(x,y,z) space.
 

From bc2f23d1a8d3798c64753ed1e6a01505361b750a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Mon, 16 Sep 2013 07:59:57 +0200
Subject: [PATCH 17/29] Revert default value of maximum iterations

---
 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 4c4a7921..2b2d8441 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -10,7 +10,7 @@ from skimage.segmentation._slic import _slic_cython
 from skimage.color import rgb2lab
 
 
-def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=0,
+def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0,
          multichannel=True, convert2lab=True, ratio=None):
     """Segments image using k-means clustering in Color-(x,y,z) space.
 

From ea1566fffb2cd08e7632be8032fd60977f30a1e2 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 16:02:24 +1000
Subject: [PATCH 18/29] Fix image dimension sanitizing at function start

`np.atleast_3d` will add a singleton dimension at the end of an array
if needed. This is not the correct thing to do if `multichannel=False`
based on the subsequent lines. If the input image was 2D with shape
`(40, 50)` and `multichannel=False`, then `np.atleast_3d` gives it
shape `(40, 50, 1)`, and then, because `multichannel=False`, the rest
of the code gives it shape `(40, 50, 1, 1)`. This results in the final
returned array having shape `(40, 50, 1)` instead of the desired
`(40, 50)`.

This commit fixes that and updates the test to detect this failure.
---
 skimage/segmentation/slic_superpixels.py | 31 ++++++++++++------------
 skimage/segmentation/tests/test_slic.py  |  2 ++
 2 files changed, 18 insertions(+), 15 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 4c4a7921..c088e17a 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -51,14 +51,12 @@ def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=0,
     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
+            - the image dimension is not 3 or 4 and `multichannel == True`
 
     Notes
     -----
-    If `sigma > 0` as is default, the image is smoothed using a Gaussian kernel
-    prior to segmentation.
+    If `sigma > 0`, the image is smoothed using a Gaussian kernel prior to
+    segmentation.
 
     The image is rescaled to be in [0, 1] prior to processing.
 
@@ -88,15 +86,18 @@ def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=0,
         compactness = ratio
 
     image = img_as_float(image)
-    image = np.atleast_3d(image)
-
-    if image.ndim == 3:
-        if multichannel:
-            # Make 2D image 3D with depth = 1
-            image = image[np.newaxis, ...]
-        else:
-            # Add channel as single last dimension
-            image = image[..., np.newaxis]
+    is2d = False
+    if image.ndim == 2:
+        # 2D grayscale image
+        image = image[np.newaxis, ..., np.newaxis]
+        is2d = True
+    elif image.ndim == 3 and multichannel:
+        # Make 2D multichannel image 3D with depth = 1
+        image = image[np.newaxis, ...]
+        is2d = True
+    elif image.ndim == 3 and not multichannel:
+        # Add channel as single last dimension
+        image = image[..., np.newaxis]
 
     if not isinstance(sigma, coll.Iterable):
         sigma = np.array([sigma, sigma, sigma])
@@ -135,7 +136,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=20, sigma=0,
 
     labels = _slic_cython(image, segments, max_iter)
 
-    if labels.shape[0] == 1:
+    if is2d:
         labels = labels[0]
 
     return labels
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index cb16e1b7..6d00716f 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -20,6 +20,7 @@ def test_color_2d():
 
     # we expect 4 segments
     assert_equal(len(np.unique(seg)), 4)
+    assert_equal(seg.shape, img.shape[:-1])
     assert_array_equal(seg[:10, :10], 0)
     assert_array_equal(seg[10:, :10], 2)
     assert_array_equal(seg[:10, 10:], 1)
@@ -39,6 +40,7 @@ def test_gray_2d():
                multichannel=False, convert2lab=False)
 
     assert_equal(len(np.unique(seg)), 4)
+    assert_equal(seg.shape, img.shape)
     assert_array_equal(seg[:10, :10], 0)
     assert_array_equal(seg[10:, :10], 2)
     assert_array_equal(seg[:10, 10:], 1)

From 88d84c720ade54dd6f3eb69a15a428fc93d1fd59 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Mon, 16 Sep 2013 08:48:46 +0200
Subject: [PATCH 19/29] Raise warning for sigma parameter and add to TODO for
 next release

---
 TODO.txt                                 | 18 +++++++++---------
 skimage/segmentation/slic_superpixels.py | 10 +++++++---
 2 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/TODO.txt b/TODO.txt
index 7ad40875..502961c8 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -1,15 +1,15 @@
 Version 0.10
 ------------
-* Remove deprecated functions:
-    - ``skimage.filter.rank.*``
-* Remove deprecated parameter ``epsilon`` of ``skimage.viewer.LineProfile``
-* Remove backwards-compatability of ``skimage.measure.regionprops``
+* Remove deprecated functions in `skimage.filter.rank.*`
+* Remove deprecated parameter `epsilon` of `skimage.viewer.LineProfile`
+* Remove backwards-compatability of `skimage.measure.regionprops`
+* Remove {`ratio`, `sigma`} deprecation warnings of `skimage.segmentation.slic`
 
 Version 0.9
 -----------
 * Remove deprecated functions
-    - ``skimage.filter.denoise_tv_chambolle``
-    - ``skimage.morphology.is_local_maximum``
-    - ``skimage.transform.hough``
-    - ``skimage.transform.probabilistic_hough``
-    - ``skimage.transform.hough_peaks``
+    - `skimage.filter.denoise_tv_chambolle`
+    - `skimage.morphology.is_local_maximum`
+    - `skimage.transform.hough`
+    - `skimage.transform.probabilistic_hough`
+    - `skimage.transform.hough_peaks`
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 010d0c62..03017cf6 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -10,7 +10,7 @@ from skimage.segmentation._slic import _slic_cython
 from skimage.color import rgb2lab
 
 
-def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0,
+def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
          multichannel=True, convert2lab=True, ratio=None):
     """Segments image using k-means clustering in Color-(x,y,z) space.
 
@@ -80,9 +80,13 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0,
     >>> segments = slic(img, n_segments=100, ratio=20)
     """
 
+    if sigma is None:
+        warnings.warn('Default value of keyword `sigma` changed from ``1`` '
+                      'to ``0``.')
+        sigma = 0
     if ratio is not None:
-        msg = 'Keyword `ratio` is deprecated. Use `compactness` instead.'
-        warnings.warn(msg)
+        warnings.warn('Keyword `ratio` is deprecated. Use `compactness` '
+                      'instead.')
         compactness = ratio
 
     image = img_as_float(image)

From bf5f08e89487983bbea771a5bc8c604066f90aaf Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 17:30:48 +1000
Subject: [PATCH 20/29] Update SLIC docstring to remove deprecated example

The `ratio` keyword has been deprecated but it was still being used in
the example in the SLIC docstring. This replaces that usage by the new
`compactness` keyword.
---
 skimage/segmentation/slic_superpixels.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 03017cf6..07bb54a5 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -75,9 +75,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     >>> 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)
+    >>> segments = slic(img, n_segments=100, compactness=10)
+    >>> # Increasing the compactness parameter yields more square regions
+    >>> segments = slic(img, n_segments=100, compactness=20)
     """
 
     if sigma is None:

From 9d86c9a1e1199c223ac435cb16a39c230c3de313 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 17:33:08 +1000
Subject: [PATCH 21/29] Ignore `convert2lab` keyword if not multichannel

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

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 07bb54a5..08ea7126 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -109,9 +109,8 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
 
-    if convert2lab:
-
-        if not multichannel or image.shape[3] != 3:
+    if convert2lab and multichannel:
+        if image.shape[3] != 3:
             raise ValueError("Lab colorspace conversion requires a RGB image.")
         image = rgb2lab(image)
 

From 610a0d1793f1b3fb91544c1921311ecf20ff355d Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 17:49:07 +1000
Subject: [PATCH 22/29] Add support for list sigma input in SLIC

Previously, having a different `sigma` for different dimensions
required an array input. This allows the user to use a simple list,
which gets converted to an array internally.

Importantly, it removes a very unhelpful error:

```python
>>> im = np.random.rand(10, 20)
>>> from skimage import segmentation as seg
Exception AttributeError: "'UmfpackContext' object has no attribute '_symbolic'" in <bound method UmfpackContext.__del__ of <scipy.sparse.linalg.dsolve.umfpack.umfpack.UmfpackContext object at 0x1045ff5d0>> ignored

>>> s = seg.slic(im, 2, sigma=[2, 1])
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-3-689b36a2f0ef> in <module>()
----> 1 s = seg.slic(im, 2, sigma=[2, 1])

/Users/nuneziglesiasj/venv/skimdev2/lib/python2.7/site-packages/scikit_image-0.9dev-py2.7-macosx-10.5-x86_64.egg/skimage/segmentation/slic_superpixels.pyc in slic(image, n_segments, compactness, max_iter, sigma, multichannel, convert2lab, ratio)
    106     if not isinstance(sigma, coll.Iterable):
    107         sigma = np.array([sigma, sigma, sigma])
--> 108     if (sigma > 0).any():
    109         sigma = list(sigma) + [0]
    110         image = ndimage.gaussian_filter(image, sigma)

AttributeError: 'bool' object has no attribute 'any'
```
---
 skimage/segmentation/slic_superpixels.py |  6 ++++--
 skimage/segmentation/tests/test_slic.py  | 11 +++++++++++
 2 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 08ea7126..be07fafe 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -27,7 +27,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         infinity, superpixel shapes become square/cubic.
     max_iter : int, optional
         Maximum number of iterations of k-means.
-    sigma : float or (3,) array of floats, optional
+    sigma : float or (3,) array-like of floats, optional
         Width of Gaussian smoothing kernel for pre-processing for each
         dimension of the image. The same sigma is applied to each dimension in
         case of a scalar value. Zero means no smoothing.
@@ -104,7 +104,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         image = image[..., np.newaxis]
 
     if not isinstance(sigma, coll.Iterable):
-        sigma = np.array([sigma, sigma, sigma])
+        sigma = np.array([sigma, sigma, sigma], float)
+    elif type(sigma) in [list, tuple]:
+        sigma = np.array(sigma, float)
     if (sigma > 0).any():
         sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 6d00716f..2a190de8 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -90,6 +90,17 @@ def test_gray_3d():
         assert_array_equal(seg[s], c)
 
 
+def test_list_sigma():
+    rnd = np.random.RandomState(0)
+    img = np.array([[1, 1, 1, 0, 0, 0],
+                    [0, 0, 0, 1, 1, 1]], np.float)
+    img += 0.1 * rnd.normal(size=img.shape)
+    result_sigma = np.array([[0, 0, 0, 1, 1, 1],
+                             [0, 0, 0, 1, 1, 1]], np.int)
+    seg_sigma = slic(img, n_segments=2, sigma=[1, 50, 1], multichannel=False)
+    assert_equal(seg_sigma, result_sigma)
+
+
 if __name__ == '__main__':
     from numpy import testing
     testing.run_module_suite()

From 846765e5f9e35cd9be24d0a391a4fde106ba6b9e Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 17:52:25 +1000
Subject: [PATCH 23/29] Add spacing support for new, speeded-up SLIC

---
 skimage/segmentation/_slic.pyx           | 15 +++++++++++----
 skimage/segmentation/slic_superpixels.py | 12 ++++++++++--
 skimage/segmentation/tests/test_slic.py  | 18 ++++++++++++++++++
 3 files changed, 39 insertions(+), 6 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index d818a501..57c9990a 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -12,7 +12,8 @@ from skimage.util import regular_grid
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
                  double[:, ::1] segments,
-                 Py_ssize_t max_iter):
+                 Py_ssize_t max_iter,
+                 double[:] spacing):
     """Helper function for SLIC segmentation.
 
     Parameters
@@ -23,6 +24,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         The initial centroids obtained by SLIC as [Z, Y, X, C...].
     max_iter : int
         The maximum number of k-means iterations.
+    spacing : 1D array of double, shape (3,)
 
     Returns
     -------
@@ -55,6 +57,11 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     cdef char change
     cdef double dist_center, cx, cy, cz, dy, dz
 
+    cdef double sz, sy, sx
+    sz = spacing[0]
+    sy = spacing[1]
+    sx = spacing[2]
+
     for i in range(max_iter):
         change = 0
         distance[:, :, :] = DBL_MAX
@@ -76,11 +83,11 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
             x_max = <Py_ssize_t>min(cx + 2 * step_x + 1, width)
 
             for z in range(z_min, z_max):
-                dz = (cz - z) ** 2
+                dz = (sz * (cz - z)) ** 2
                 for y in range(y_min, y_max):
-                    dy = (cy - y) ** 2
+                    dy = (sy * (cy - y)) ** 2
                     for x in range(x_min, x_max):
-                        dist_center = dz + dy + (cx - x) ** 2
+                        dist_center = dz + dy + (sx * (cx - x)) ** 2
                         for c in range(3, n_features):
                             dist_center += (image_zyx[z, y, x, c - 3]
                                             - segments[k, c]) ** 2
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index be07fafe..e9ae8ee7 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -11,7 +11,7 @@ from skimage.color import rgb2lab
 
 
 def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
-         multichannel=True, convert2lab=True, ratio=None):
+         spacing=None, multichannel=True, convert2lab=True, ratio=None):
     """Segments image using k-means clustering in Color-(x,y,z) space.
 
     Parameters
@@ -31,6 +31,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         Width of Gaussian smoothing kernel for pre-processing for each
         dimension of the image. The same sigma is applied to each dimension in
         case of a scalar value. Zero means no smoothing.
+    spacing : (3,) array-like of floats, optional
+        The voxel spacing along each image dimension. By default, `slic`
+        assumes uniform spacing (same voxel resolution along z, y and x).
     multichannel : bool, optional
         Whether the last axis of the image is to be interpreted as multiple
         channels or another spatial dimension.
@@ -103,11 +106,16 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         # Add channel as single last dimension
         image = image[..., np.newaxis]
 
+    if spacing is None:
+        spacing = np.ones(3)
+    elif type(spacing) in [list, tuple]:
+        spacing = np.array(spacing, float)
     if not isinstance(sigma, coll.Iterable):
         sigma = np.array([sigma, sigma, sigma], float)
     elif type(sigma) in [list, tuple]:
         sigma = np.array(sigma, float)
     if (sigma > 0).any():
+        sigma /= spacing.astype(float)
         sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
 
@@ -139,7 +147,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     ratio = float(max((step_z, step_y, step_x))) / compactness
     image = np.ascontiguousarray(image * ratio)
 
-    labels = _slic_cython(image, segments, max_iter)
+    labels = _slic_cython(image, segments, max_iter, spacing)
 
     if is2d:
         labels = labels[0]
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 2a190de8..a4657785 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -101,6 +101,24 @@ def test_list_sigma():
     assert_equal(seg_sigma, result_sigma)
 
 
+def test_spacing():
+    rnd = np.random.RandomState(0)
+    img = np.array([[1, 1, 1, 0, 0],
+                    [1, 1, 0, 0, 0]], np.float)
+    result_non_spaced = np.array([[0, 0, 0, 1, 1],
+                                  [0, 0, 1, 1, 1]], np.int)
+    result_spaced = np.array([[0, 0, 0, 0, 0],
+                              [1, 1, 1, 1, 1]], np.int)
+    img += 0.1 * rnd.normal(size=img.shape)
+    seg_non_spaced = slic(img, n_segments=2, sigma=0, multichannel=False,
+                          compactness=1.0)
+    seg_spaced = slic(img, n_segments=2, sigma=0, spacing=[1, 500, 1],
+                      compactness=1.0, multichannel=False)
+    assert_equal(seg_non_spaced, result_non_spaced)
+    assert_equal(seg_spaced, result_spaced)
+
+
+
 if __name__ == '__main__':
     from numpy import testing
     testing.run_module_suite()

From 00e5ff263bc1f32b41f02f9ca8dbafbb05459668 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 21:26:44 +1000
Subject: [PATCH 24/29] Add `spacing` descriptions; use np.double

Implement @ahojnnes's comments on pull request. Use `np.double` as
dtype for arrays because in Cython, `float` is not `np.double`. And
add further clarification about the `spacing` parameter.
---
 skimage/segmentation/_slic.pyx           |  3 +++
 skimage/segmentation/slic_superpixels.py | 19 +++++++++++++------
 2 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 57c9990a..3247dd79 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -25,6 +25,9 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     max_iter : int
         The maximum number of k-means iterations.
     spacing : 1D array of double, shape (3,)
+        The voxel spacing along each image dimension. This parameter
+        controls the weights of the distances along z, y, and x during
+        k-means clustering.
 
     Returns
     -------
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index e9ae8ee7..3cb67874 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -34,6 +34,8 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     spacing : (3,) array-like of floats, optional
         The voxel spacing along each image dimension. By default, `slic`
         assumes uniform spacing (same voxel resolution along z, y and x).
+        This parameter controls the weights of the distances along z, y,
+        and x during k-means clustering.
     multichannel : bool, optional
         Whether the last axis of the image is to be interpreted as multiple
         channels or another spatial dimension.
@@ -61,6 +63,11 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     If `sigma > 0`, the image is smoothed using a Gaussian kernel prior to
     segmentation.
 
+    If `sigma > 0` and `spacing` is provided, the kernel width is divided
+    along each dimension by the spacing. For example, if `sigma=1` and
+    `spacing=[5, 1, 1]`, the effective `sigma` is `[0.2, 1, 1]`. This
+    ensures sensible smoothing for anisotropic images.
+
     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
@@ -108,14 +115,14 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
 
     if spacing is None:
         spacing = np.ones(3)
-    elif type(spacing) in [list, tuple]:
-        spacing = np.array(spacing, float)
+    elif isinstance(spacing, (list, tuple)):
+        spacing = np.array(spacing, np.double)
     if not isinstance(sigma, coll.Iterable):
-        sigma = np.array([sigma, sigma, sigma], float)
-    elif type(sigma) in [list, tuple]:
-        sigma = np.array(sigma, float)
+        sigma = np.array([sigma, sigma, sigma], np.double)
+    elif isinstance(spacing, (list, tuple)):
+        sigma = np.array(sigma, np.double)
     if (sigma > 0).any():
-        sigma /= spacing.astype(float)
+        sigma /= spacing.astype(np.double)
         sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
 

From 05fbc3fbfcc00157841c18b15b171d6477bfb5da Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Mon, 16 Sep 2013 22:42:48 +1000
Subject: [PATCH 25/29] Bug fix: typo: wrote spacing instead of sigma

---
 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 3cb67874..422dff98 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -119,7 +119,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         spacing = np.array(spacing, np.double)
     if not isinstance(sigma, coll.Iterable):
         sigma = np.array([sigma, sigma, sigma], np.double)
-    elif isinstance(spacing, (list, tuple)):
+    elif isinstance(sigma, (list, tuple)):
         sigma = np.array(sigma, np.double)
     if (sigma > 0).any():
         sigma /= spacing.astype(np.double)

From a086f03f39187858bbd65389a2a4b2e302d37e7c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 19 Sep 2013 08:36:56 +0200
Subject: [PATCH 26/29] Document default value API change

---
 doc/source/api_changes.txt | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/doc/source/api_changes.txt b/doc/source/api_changes.txt
index c363288d..bf9ddb60 100644
--- a/doc/source/api_changes.txt
+++ b/doc/source/api_changes.txt
@@ -1,6 +1,8 @@
 Version 0.9
 -----------
-- No longer wrap ``imread`` output in an ``Image`` class.
+- No longer wrap ``imread`` output in an ``Image`` class
+- Change default value of `sigma` parameter in ``skimage.segmentation.slic``
+  to 0
 
 Version 0.4
 -----------

From e17e7699b328c2489b4aea0a8aacd312b5e90a12 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Sat, 21 Sep 2013 00:12:12 +1000
Subject: [PATCH 27/29] Add note about (z, y, x) order in slic cython doc

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 3247dd79..36cfa1a5 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -33,6 +33,27 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     -------
     nearest_segments : 3D array of int, shape (Z, Y, X)
         The label field/superpixels found by SLIC.
+
+    Notes
+    -----
+    The image is considered to be in (z, y, x) order, which can be
+    surprising. More commonly, the order (x, y, z) is used. However,
+    in 3D image analysis, 'z' is usually the "special" dimension, with,
+    for example, a different effective resolution than the other two
+    axes. Therefore, x and y are often processed together, or viewed as
+    a cut-plane through the volume. So, if the order was (x, y, z) and
+    we wanted to look at the 5th cut plane, we would write::
+
+        my_z_plane = img3d[:, :, 5]
+
+    but, assuming a C-contiguous array, this would grab a discontiguous
+    slice of memory, which is bad for performance. In contrast, if we
+    see the image as (z, y, x) ordered, we would do::
+
+        my_z_plane = img3d[5]
+
+    and get back a contiguous block of memory. This is better both for
+    performance and for readability.
     """
 
     # initialize on grid

From 845d083470c065cfb1316e788ebb160abd222e12 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sun, 22 Sep 2013 07:22:09 +0200
Subject: [PATCH 28/29] Define 1d arrays as contiguous

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

diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 36cfa1a5..c3d95ee0 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -13,7 +13,7 @@ from skimage.util import regular_grid
 def _slic_cython(double[:, :, :, ::1] image_zyx,
                  double[:, ::1] segments,
                  Py_ssize_t max_iter,
-                 double[:] spacing):
+                 double[::1] spacing):
     """Helper function for SLIC segmentation.
 
     Parameters
@@ -75,7 +75,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         = np.empty((depth, height, width), dtype=np.intp)
     cdef double[:, :, ::1] distance \
         = np.empty((depth, height, width), dtype=np.double)
-    cdef Py_ssize_t[:] n_segment_elems = np.zeros(n_segments, dtype=np.intp)
+    cdef Py_ssize_t[::1] n_segment_elems = np.zeros(n_segments, dtype=np.intp)
 
     cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max
     cdef char change

From 8fb6fd02a0975941f888c9271e49f65483298f29 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 26 Sep 2013 16:49:55 +0200
Subject: [PATCH 29/29] Scale sigma only in scalar case

---
 skimage/segmentation/slic_superpixels.py | 42 +++++++++++++-----------
 1 file changed, 23 insertions(+), 19 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 422dff98..8276b1a8 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -31,6 +31,8 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         Width of Gaussian smoothing kernel for pre-processing for each
         dimension of the image. The same sigma is applied to each dimension in
         case of a scalar value. Zero means no smoothing.
+        Note, that `sigma` is automatically scaled if it is scalar and a
+        manual voxel spacing is provided (see Notes section).
     spacing : (3,) array-like of floats, optional
         The voxel spacing along each image dimension. By default, `slic`
         assumes uniform spacing (same voxel resolution along z, y and x).
@@ -60,19 +62,19 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
 
     Notes
     -----
-    If `sigma > 0`, the image is smoothed using a Gaussian kernel prior to
-    segmentation.
+    * If `sigma > 0`, the image is smoothed using a Gaussian kernel prior to
+      segmentation.
 
-    If `sigma > 0` and `spacing` is provided, the kernel width is divided
-    along each dimension by the spacing. For example, if `sigma=1` and
-    `spacing=[5, 1, 1]`, the effective `sigma` is `[0.2, 1, 1]`. This
-    ensures sensible smoothing for anisotropic images.
+    * If `sigma` is scalar and `spacing` is provided, the kernel width is
+      divided along each dimension by the spacing. For example, if ``sigma=1``
+      and ``spacing=[5, 1, 1]``, the effective `sigma` is ``[0.2, 1, 1]``. This
+      ensures sensible smoothing for anisotropic images.
 
-    The image is rescaled to be in [0, 1] prior to processing.
+    * 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`.
+    * 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
     ----------
@@ -100,15 +102,15 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         compactness = ratio
 
     image = img_as_float(image)
-    is2d = False
+    is_2d = False
     if image.ndim == 2:
         # 2D grayscale image
         image = image[np.newaxis, ..., np.newaxis]
-        is2d = True
+        is_2d = True
     elif image.ndim == 3 and multichannel:
         # Make 2D multichannel image 3D with depth = 1
         image = image[np.newaxis, ...]
-        is2d = True
+        is_2d = True
     elif image.ndim == 3 and not multichannel:
         # Add channel as single last dimension
         image = image[..., np.newaxis]
@@ -116,13 +118,15 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     if spacing is None:
         spacing = np.ones(3)
     elif isinstance(spacing, (list, tuple)):
-        spacing = np.array(spacing, np.double)
+        spacing = np.array(spacing, dtype=np.double)
+
     if not isinstance(sigma, coll.Iterable):
-        sigma = np.array([sigma, sigma, sigma], np.double)
-    elif isinstance(sigma, (list, tuple)):
-        sigma = np.array(sigma, np.double)
-    if (sigma > 0).any():
+        sigma = np.array([sigma, sigma, sigma], dtype=np.double)
         sigma /= spacing.astype(np.double)
+    elif isinstance(sigma, (list, tuple)):
+        sigma = np.array(sigma, dtype=np.double)
+    if (sigma > 0).any():
+        # add zero smoothing for multichannel dimension
         sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
 
@@ -156,7 +160,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
 
     labels = _slic_cython(image, segments, max_iter, spacing)
 
-    if is2d:
+    if is_2d:
         labels = labels[0]
 
     return labels