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)