diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx
index 64d09520..f7caee03 100644
--- a/skimage/segmentation/_slic.pyx
+++ b/skimage/segmentation/_slic.pyx
@@ -3,16 +3,20 @@
 #cython: nonecheck=False
 #cython: wraparound=False
 from libc.float cimport DBL_MAX
+from cpython cimport bool
 
 import numpy as np
 cimport numpy as cnp
 
 from skimage.util import regular_grid
 
+
 def _slic_cython(double[:, :, :, ::1] image_zyx,
                  double[:, ::1] segments,
+                 float step,
                  Py_ssize_t max_iter,
-                 double[::1] spacing):
+                 double[::1] spacing,
+                 bint slic_zero):
     """Helper function for SLIC segmentation.
 
     Parameters
@@ -21,12 +25,16 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         The input image.
     segments : 2D array of double, shape (N, 3 + C)
         The initial centroids obtained by SLIC as [Z, Y, X, C...].
+    step : double
+        The size of the step between two seeds in voxels.
     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.
+        k-means clustering.        
+    slic_zero : bool
+        True to run SLIC-zero, False to run original SLIC.
 
     Returns
     -------
@@ -85,6 +93,14 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     sy = spacing[1]
     sx = spacing[2]
 
+    # The colors are scaled before being passed to _slic_cython so
+    # max_color_sq can be initialised as all ones
+    cdef double[::1] max_dist_color = np.ones(n_segments, dtype=np.double)
+    cdef double dist_color
+
+    # The reference implementation (Achanta et al.) calls this invxywt
+    cdef double spatial_weight = float(1) / (step ** 2)
+
     for i in range(max_iter):
         change = 0
         distance[:, :, :] = DBL_MAX
@@ -110,10 +126,16 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
                 for y in range(y_min, y_max):
                     dy = (sy * (cy - y)) ** 2
                     for x in range(x_min, x_max):
-                        dist_center = dz + dy + (sx * (cx - x)) ** 2
+                        dist_center = (dz + dy + (sx * (cx - x)) ** 2) * spatial_weight
+                        dist_color = 0
                         for c in range(3, n_features):
-                            dist_center += (image_zyx[z, y, x, c - 3]
+                            dist_color += (image_zyx[z, y, x, c - 3]
                                             - segments[k, c]) ** 2
+                        if slic_zero:
+                            dist_center += dist_color / max_dist_color[k]
+                        else:
+                            dist_center += dist_color
+
                         if distance[z, y, x] > dist_center:
                             nearest_segments[z, y, x] = k
                             distance[z, y, x] = dist_center
@@ -144,6 +166,24 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
             for c in range(n_features):
                 segments[k, c] /= n_segment_elems[k]
 
+        # If in SLICO mode, update the color distance maxima
+        if slic_zero:
+            for z in range(depth):
+                for y in range(height):
+                    for x in range(width):
+
+                        k = nearest_segments[z, y, x]
+                        dist_color = 0
+
+                        for c in range(3, n_features):
+                            dist_color += (image_zyx[z, y, x, c - 3] -
+                                           segments[k, c]) ** 2
+                        
+                        # The reference implementation seems to only change
+                        # the color if it increases from previous iteration
+                        if max_dist_color[k] < dist_color:
+                            max_dist_color[k] = dist_color
+
     return np.asarray(nearest_segments)
 
 
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 9be5aa62..11373381 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -12,7 +12,8 @@ from skimage.color import rgb2lab
 
 def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
          spacing=None, multichannel=True, convert2lab=True, ratio=None,
-         enforce_connectivity=False, min_size_factor=0.5, max_size_factor=3):
+         enforce_connectivity=False, min_size_factor=0.5, max_size_factor=3,
+         slic_zero=False):
     """Segments image using k-means clustering in Color-(x,y,z) space.
 
     Parameters
@@ -25,7 +26,8 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     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.
+        infinity, superpixel shapes become square/cubic. In SLICO mode, this
+        is the initial compactness.
     max_iter : int, optional
         Maximum number of iterations of k-means.
     sigma : float or (3,) array-like of floats, optional
@@ -56,6 +58,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     max_size_factor: float, optional
         Proportion of the maximum connected segment size. A value of 3 works
         in most of the cases.
+    slic_zero: bool, optional
+        Run SLIC-zero, the zero-parameter mode of SLIC
+
     Returns
     -------
     labels : 2D or 3D array
@@ -165,16 +170,18 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     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])
+                               segments_color],
+                              axis=-1).reshape(-1, 3 + image.shape[3])
     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
+    step = float(max((step_z, step_y, step_x)))
+    ratio = 1.0 / compactness
+
     image = np.ascontiguousarray(image * ratio)
 
-    labels = _slic_cython(image, segments, max_iter, spacing)
+    labels = _slic_cython(image, segments, step, max_iter, spacing, slic_zero)
 
     if enforce_connectivity:
         segment_size = depth * height * width / n_segments
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index 15f03437..7dda66d2 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -147,6 +147,30 @@ def test_enforce_connectivity():
     assert_equal(segments_connected, result_connected)
     assert_equal(segments_disconnected, result_disconnected)
 
+
+def test_slic_zero():
+    # Same as test_color_2d but with slic_zero=True
+    rnd = np.random.RandomState(0)
+    img = np.zeros((20, 21, 3))
+    img[:10, :10, 0] = 1
+    img[10:, :10, 1] = 1
+    img[10:, 10:, 2] = 1
+    img += 0.01 * rnd.normal(size=img.shape)
+    img[img > 1] = 1
+    img[img < 0] = 0
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        seg = slic(img, n_segments=4, sigma=0, slic_zero=True)
+
+    # we expect 4 segments
+    assert_equal(len(np.unique(seg)), 4)
+    assert_equal(seg.shape, img.shape[:-1])
+    assert_equal(seg[:10, :10], 0)
+    assert_equal(seg[10:, :10], 2)
+    assert_equal(seg[:10, 10:], 1)
+    assert_equal(seg[10:, 10:], 3)
+
+
 if __name__ == '__main__':
     from numpy import testing