diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index a4274484..474b2966 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -11,6 +11,7 @@ if chk < 18:    # Use internal version for numpy versions < 1.8.x
 else:
     from numpy import pad
 del numpy, ver, chk
+from .regular_grid import regular_grid
 
 
 __all__ = ['img_as_float',
@@ -23,3 +24,4 @@ __all__ = ['img_as_float',
            'view_as_windows',
            'pad',
            'random_noise']
+           'regular_grid']
diff --git a/skimage/util/regular_grid.py b/skimage/util/regular_grid.py
new file mode 100644
index 00000000..40980491
--- /dev/null
+++ b/skimage/util/regular_grid.py
@@ -0,0 +1,44 @@
+import numpy as np
+
+def regular_grid(ar_shape, n_points):
+    """Find `n_points` regularly spaced along `ar_shape`.
+    
+    The returned points (as slices) should be as close to cubically-spaced as
+    possible.
+
+    Parameters
+    ----------
+    ar_shape : array-like of ints
+        The shape of the space embedding the grid. `len(ar_shape)` is the
+        number of dimensions.
+    n_points : int
+        The (approximate) number of points to embed in the space.
+
+    Returns
+    -------
+    slices : list of slice objects
+        A slice along each dimension of `ar_shape`, such that the intersection
+        of all the slices give the coordinates of regularly spaced points.
+    """
+    ar_shape = np.asanyarray(ar_shape)
+    ndim = len(ar_shape)
+    unsort_dim_idxs = np.argsort(np.argsort(ar_shape))
+    sorted_dims = np.sort(ar_shape)
+    space_size = float(np.prod(ar_shape))
+    if space_size <= n_points:
+        return [slice(None)] * ndim
+    stepsizes = (space_size / n_points) ** (1.0 / ndim) * np.ones(ndim)
+    if (sorted_dims < stepsizes).any():
+        for dim in range(ndim):
+            stepsizes[dim] = sorted_dims[dim]
+            space_size = float(np.prod(sorted_dims[dim+1:]))
+            stepsizes[dim+1:] = ((space_size / n_points) **
+                                            (1.0 / (ndim - dim - 1)))
+            if (sorted_dims >= stepsizes).all():
+                break
+    starts = np.floor(stepsizes/2)
+    stepsizes = np.round(stepsizes)
+    slices = [slice(start, None, step) for 
+                                        start, step in zip(starts, stepsizes)]
+    slices = [slices[i] for i in unsort_dim_idxs]
+    return slices