diff --git a/bento.info b/bento.info
index 6a8a0000..a36a4ebb 100644
--- a/bento.info
+++ b/bento.info
@@ -79,9 +79,6 @@ Library:
     Extension: skimage.graph._spath
         Sources:
             skimage/graph/_spath.pyx
-    Extension: skimage.morphology.cmorph
-        Sources:
-            skimage/morphology/cmorph.pyx
     Extension: skimage.graph.heap
         Sources:
             skimage/graph/heap.pyx
diff --git a/skimage/filters/rank/README.rst b/skimage/filters/rank/README.rst
index e5c5a9ad..6e063ec2 100644
--- a/skimage/filters/rank/README.rst
+++ b/skimage/filters/rank/README.rst
@@ -16,7 +16,7 @@ followed by the moving window is given hereunder
 /--------------------------/
 \-------------------------- ...
 
-We compare cmorph.dilate to this histogram based method to show how
+We compare grey.dilate to this histogram based method to show how
 computational costs increase with respect to image size or structuring element
 size. This implementation gives better results for large structuring elements.
 
@@ -26,7 +26,7 @@ update the local histogram. The histogram size is 8-bit (256 bins) for 8-bit
 images and 2 to 16-bit for 16-bit images depending on the maximum value of the
 image.
 
-The filter is applied up to the image border, the neighboorhood used is
+The filter is applied up to the image border, the neighborhood used is
 adjusted accordingly. The user may provide a mask image (same size as input
 image) where non zero values are the part of the image participating in the
 histogram computation. By default the entire image is filtered.
diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py
index 64c3ca97..48628432 100644
--- a/skimage/filters/rank/tests/test_rank.py
+++ b/skimage/filters/rank/tests/test_rank.py
@@ -5,7 +5,7 @@ from numpy.testing import run_module_suite, assert_equal, assert_raises
 import skimage
 from skimage import img_as_ubyte, img_as_float
 from skimage import data, util, morphology
-from skimage.morphology import cmorph, disk
+from skimage.morphology import grey, disk
 from skimage.filters import rank
 from skimage._shared._warnings import expected_warnings
 
@@ -87,7 +87,6 @@ def check_all():
 def test_random_sizes():
     # make sure the size is not a problem
 
-    niter = 10
     elem = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype=np.uint8)
     for m, n in np.random.random_integers(1, 100, size=(10, 2)):
         mask = np.ones((m, n), dtype=np.uint8)
@@ -118,31 +117,31 @@ def test_random_sizes():
         assert_equal(image16.shape, out16.shape)
 
 
-def test_compare_with_cmorph_dilate():
+def test_compare_with_grey_dilation():
     # compare the result of maximum filter with dilate
 
     image = (np.random.rand(100, 100) * 256).astype(np.uint8)
     out = np.empty_like(image)
     mask = np.ones(image.shape, dtype=np.uint8)
 
-    for r in range(1, 20, 1):
+    for r in range(3, 20, 2):
         elem = np.ones((r, r), dtype=np.uint8)
         rank.maximum(image=image, selem=elem, out=out, mask=mask)
-        cm = cmorph._dilate(image=image, selem=elem)
+        cm = grey.dilation(image=image, selem=elem)
         assert_equal(out, cm)
 
 
-def test_compare_with_cmorph_erode():
+def test_compare_with_grey_erosion():
     # compare the result of maximum filter with erode
 
     image = (np.random.rand(100, 100) * 256).astype(np.uint8)
     out = np.empty_like(image)
     mask = np.ones(image.shape, dtype=np.uint8)
 
-    for r in range(1, 20, 1):
+    for r in range(3, 20, 2):
         elem = np.ones((r, r), dtype=np.uint8)
         rank.minimum(image=image, selem=elem, out=out, mask=mask)
-        cm = cmorph._erode(image=image, selem=elem)
+        cm = grey.erosion(image=image, selem=elem)
         assert_equal(out, cm)
 
 
diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py
index 4244738f..0ff5c8f7 100644
--- a/skimage/morphology/binary.py
+++ b/skimage/morphology/binary.py
@@ -2,15 +2,13 @@
 Binary morphological operations
 """
 import numpy as np
-from scipy import ndimage
-from .misc import default_fallback
+from scipy import ndimage as nd
+from .misc import default_selem
 
 
-# Our functions only work in 2D, so for 3D or higher input we should fall back
-# on `scipy.ndimage`. Additionally, we want to use a cross-shaped structuring
-# element of the appropriate dimension for each of these functions.
-# The `default_callback` provides all these.
-@default_fallback
+# The default_selem decorator provides a diamond structuring element as default
+# with the same dimension as the input image and size 3 along each axis.
+@default_selem
 def binary_erosion(image, selem=None, out=None):
     """Return fast binary morphological erosion of an image.
 
@@ -35,27 +33,17 @@ def binary_erosion(image, selem=None, out=None):
     Returns
     -------
     eroded : ndarray of bool or uint
-        The result of the morphological erosion with values in ``[0, 1]``.
+        The result of the morphological erosion taking values in
+        ``[False, True]``.
 
     """
-
-    selem = (selem != 0)
-    selem_sum = np.sum(selem)
-
-    if selem_sum <= 255:
-        conv = np.empty_like(image, dtype=np.uint8)
-    else:
-        conv = np.empty_like(image, dtype=np.uint)
-
-    binary = (image > 0).view(np.uint8)
-    ndimage.convolve(binary, selem, mode='constant', cval=1, output=conv)
-
     if out is None:
-        out = np.empty_like(conv, dtype=np.bool)
-    return np.equal(conv, selem_sum, out=out)
+        out = np.empty(image.shape, dtype=np.bool)
+    nd.binary_erosion(image, structure=selem, output=out)
+    return out
 
 
-@default_fallback
+@default_selem
 def binary_dilation(image, selem=None, out=None):
     """Return fast binary morphological dilation of an image.
 
@@ -81,26 +69,16 @@ def binary_dilation(image, selem=None, out=None):
     Returns
     -------
     dilated : ndarray of bool or uint
-        The result of the morphological dilation with values in ``[0, 1]``.
-
+        The result of the morphological dilation with values in
+        ``[False, True]``.
     """
-
-    selem = (selem != 0)
-
-    if np.sum(selem) <= 255:
-        conv = np.empty_like(image, dtype=np.uint8)
-    else:
-        conv = np.empty_like(image, dtype=np.uint)
-
-    binary = (image > 0).view(np.uint8)
-    ndimage.convolve(binary, selem, mode='constant', cval=0, output=conv)
-
     if out is None:
-        out = np.empty_like(conv, dtype=np.bool)
-    return np.not_equal(conv, 0, out=out)
+        out = np.empty(image.shape, dtype=np.bool)
+    nd.binary_dilation(image, structure=selem, output=out)
+    return out
 
 
-@default_fallback
+@default_selem
 def binary_opening(image, selem=None, out=None):
     """Return fast binary morphological opening of an image.
 
@@ -134,7 +112,7 @@ def binary_opening(image, selem=None, out=None):
     return out
 
 
-@default_fallback
+@default_selem
 def binary_closing(image, selem=None, out=None):
     """Return fast binary morphological closing of an image.
 
@@ -163,7 +141,6 @@ def binary_closing(image, selem=None, out=None):
         The result of the morphological closing.
 
     """
-
     dilated = binary_dilation(image, selem)
     out = binary_erosion(dilated, selem, out=out)
     return out
diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py
index e65e756b..8ab162cf 100644
--- a/skimage/morphology/grey.py
+++ b/skimage/morphology/grey.py
@@ -1,17 +1,134 @@
 """
 Grayscale morphological operations
 """
-from skimage import img_as_ubyte
-from .misc import default_fallback
-
-from . import cmorph
+import functools
 
+import numpy as np
+from scipy import ndimage as nd
+from .misc import default_selem
+from ..util import pad, crop
 
 __all__ = ['erosion', 'dilation', 'opening', 'closing', 'white_tophat',
            'black_tophat']
 
 
-@default_fallback
+def _shift_selem(selem, shift_x, shift_y):
+    """Shift the binary image `selem` in the left and/or up.
+
+    This only affects 2D structuring elements with even number of rows
+    or columns.
+
+    Parameters
+    ----------
+    selem : 2D array, shape (M, N)
+        The input structuring element.
+    shift_x, shift_y : bool
+        Whether to move `selem` along each axis.
+
+    Returns
+    -------
+    out : 2D array, shape (M + int(shift_x), N + int(shift_y))
+        The shifted structuring element.
+    """
+    if selem.ndim > 2:
+        # do nothing for 3D or higher structuring elements
+        return selem
+    m, n = selem.shape
+    if m % 2 == 0:
+        extra_row = np.zeros((1, n), selem.dtype)
+        if shift_x:
+            selem = np.vstack((selem, extra_row))
+        else:
+            selem = np.vstack((extra_row, selem))
+        m += 1
+    if n % 2 == 0:
+        extra_col = np.zeros((m, 1), selem.dtype)
+        if shift_y:
+            selem = np.hstack((selem, extra_col))
+        else:
+            selem = np.hstack((extra_col, selem))
+    return selem
+
+
+def _invert_selem(selem):
+    """Change the order of the values in `selem`.
+
+    This is a patch for the *weird* footprint inversion in
+    `nd.grey_morphology` [1].
+
+    Parameters
+    ----------
+    selem : array
+        The input structuring element.
+
+    Returns
+    -------
+    inverted : array, same shape and type as `selem`
+        The structuring element, in opposite order.
+
+    Examples
+    --------
+    >>> selem = np.array([[0, 0, 0], [0, 1, 1], [0, 1, 1]], np.uint8)
+    >>> _invert_selem(selem)
+    array([[1, 1, 0],
+           [1, 1, 0],
+           [0, 0, 0]], dtype=uint8)
+
+    References
+    ----------
+    [1] https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285
+    """
+    inverted = selem[(slice(None, None, -1),) * selem.ndim]
+    return inverted
+
+
+def pad_for_eccentric_selems(func):
+    """Pad input images for certain morphological operations.
+
+    Parameters
+    ----------
+    func : callable
+        A morphological function, either opening or closing, that
+        supports eccentric structuring elements. Its parameters must
+        include at least `image`, `selem`, and `out`.
+
+    Returns
+    -------
+    func_out : callable
+        The same function, but correctly padding the input image before
+        applying the input function.
+
+    See Also
+    --------
+    opening, closing.
+    """
+    @functools.wraps(func)
+    def func_out(image, selem, out=None, *args, **kwargs):
+        pad_widths = []
+        padding = False
+        if out is None:
+            out = np.empty_like(image)
+        for axis_len in selem.shape:
+            if axis_len % 2 == 0:
+                axis_pad_width = axis_len - 1
+                padding = True
+            else:
+                axis_pad_width = 0
+            pad_widths.append((axis_pad_width,) * 2)
+        if padding:
+            image = pad(image, pad_widths, mode='edge')
+            out_temp = np.empty_like(image)
+        else:
+            out_temp = out
+        out_temp = func(image, selem, out=out_temp, *args, **kwargs)
+        if padding:
+            out[:] = crop(out_temp, pad_widths)
+        else:
+            out = out_temp
+        return out
+    return func_out
+
+@default_selem
 def erosion(image, selem=None, out=None, shift_x=False, shift_y=False):
     """Return greyscale morphological erosion of an image.
 
@@ -24,7 +141,7 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False):
     image : ndarray
         Image array.
     selem : ndarray, optional
-        The neighborhood expressed as a 2-D array of 1's and 0's.
+        The neighborhood expressed as an array of 1's and 0's.
         If None, use cross-shaped structuring element (connectivity=1).
     out : ndarrays, optional
         The array to store the result of the morphology. If None is
@@ -35,14 +152,14 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False):
 
     Returns
     -------
-    eroded : uint8 array
+    eroded : array, same shape as `image`
         The result of the morphological erosion.
 
     Notes
     -----
-    For `uint8` (and `uint16` up to a certain bit-depth) data, the lower
-    algorithm complexity makes the `skimage.filter.rank.minimum` function more
-    efficient for larger images and structuring elements.
+    For ``uint8`` (and ``uint16`` up to a certain bit-depth) data, the
+    lower algorithm complexity makes the `skimage.filter.rank.minimum`
+    function more efficient for larger images and structuring elements.
 
     Examples
     --------
@@ -62,16 +179,15 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False):
            [0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-
-    if image is out:
-        raise NotImplementedError("In-place erosion not supported!")
-    image = img_as_ubyte(image)
-    selem = img_as_ubyte(selem)
-    return cmorph._erode(image, selem, out=out,
-                         shift_x=shift_x, shift_y=shift_y)
+    selem = np.array(selem)
+    selem = _shift_selem(selem, shift_x, shift_y)
+    if out is None:
+        out = np.empty_like(image)
+    nd.grey_erosion(image, footprint=selem, output=out)
+    return out
 
 
-@default_fallback
+@default_selem
 def dilation(image, selem=None, out=None, shift_x=False, shift_y=False):
     """Return greyscale morphological dilation of an image.
 
@@ -96,7 +212,7 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False):
 
     Returns
     -------
-    dilated : uint8 array
+    dilated : uint8 array, same shape and type as `image`
         The result of the morphological dilation.
 
     Notes
@@ -123,17 +239,22 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False):
            [0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-
-    if image is out:
-        raise NotImplementedError("In-place dilation not supported!")
-
-    image = img_as_ubyte(image)
-    selem = img_as_ubyte(selem)
-    return cmorph._dilate(image, selem, out=out,
-                          shift_x=shift_x, shift_y=shift_y)
+    selem = np.array(selem)
+    selem = _shift_selem(selem, shift_x, shift_y)
+    # Inside ndimage.grey_dilation, the structuring element is inverted,
+    # eg. `selem = selem[::-1, ::-1]` for 2D [1]_, for reasons unknown to
+    # this author (@jni). To "patch" this behaviour, we invert our own
+    # selem before passing it to `nd.grey_dilation`.
+    # [1] https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285
+    selem = _invert_selem(selem)
+    if out is None:
+        out = np.empty_like(image)
+    nd.grey_dilation(image, footprint=selem, output=out)
+    return out
 
 
-@default_fallback
+@default_selem
+@pad_for_eccentric_selems
 def opening(image, selem=None, out=None):
     """Return greyscale morphological opening of an image.
 
@@ -147,7 +268,7 @@ def opening(image, selem=None, out=None):
     image : ndarray
         Image array.
     selem : ndarray, optional
-        The neighborhood expressed as a 2-D array of 1's and 0's.
+        The neighborhood expressed as an array of 1's and 0's.
         If None, use cross-shaped structuring element (connectivity=1).
     out : ndarray, optional
         The array to store the result of the morphology. If None
@@ -155,7 +276,7 @@ def opening(image, selem=None, out=None):
 
     Returns
     -------
-    opening : uint8 array
+    opening : array, same shape and type as `image`
         The result of the morphological opening.
 
     Examples
@@ -176,17 +297,14 @@ def opening(image, selem=None, out=None):
            [0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-
-    h, w = selem.shape
-    shift_x = True if (w % 2) == 0 else False
-    shift_y = True if (h % 2) == 0 else False
-
     eroded = erosion(image, selem)
-    out = dilation(eroded, selem, out=out, shift_x=shift_x, shift_y=shift_y)
+    # note: shift_x, shift_y do nothing if selem side length is odd
+    out = dilation(eroded, selem, out=out, shift_x=True, shift_y=True)
     return out
 
 
-@default_fallback
+@default_selem
+@pad_for_eccentric_selems
 def closing(image, selem=None, out=None):
     """Return greyscale morphological closing of an image.
 
@@ -200,7 +318,7 @@ def closing(image, selem=None, out=None):
     image : ndarray
         Image array.
     selem : ndarray, optional
-        The neighborhood expressed as a 2-D array of 1's and 0's.
+        The neighborhood expressed as an array of 1's and 0's.
         If None, use cross-shaped structuring element (connectivity=1).
     out : ndarray, optional
         The array to store the result of the morphology. If None,
@@ -208,7 +326,7 @@ def closing(image, selem=None, out=None):
 
     Returns
     -------
-    closing : uint8 array
+    closing : array, same shape and type as `image`
         The result of the morphological closing.
 
     Examples
@@ -229,17 +347,13 @@ def closing(image, selem=None, out=None):
            [0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-
-    h, w = selem.shape
-    shift_x = True if (w % 2) == 0 else False
-    shift_y = True if (h % 2) == 0 else False
-
     dilated = dilation(image, selem)
-    out = erosion(dilated, selem, out=out, shift_x=shift_x, shift_y=shift_y)
+    # note: shift_x, shift_y do nothing if selem side length is odd
+    out = erosion(dilated, selem, out=out, shift_x=True, shift_y=True)
     return out
 
 
-@default_fallback
+@default_selem
 def white_tophat(image, selem=None, out=None):
     """Return white top hat of an image.
 
@@ -252,7 +366,7 @@ def white_tophat(image, selem=None, out=None):
     image : ndarray
         Image array.
     selem : ndarray, optional
-        The neighborhood expressed as a 2-D array of 1's and 0's.
+        The neighborhood expressed as an array of 1's and 0's.
         If None, use cross-shaped structuring element (connectivity=1).
     out : ndarray, optional
         The array to store the result of the morphology. If None
@@ -260,7 +374,7 @@ def white_tophat(image, selem=None, out=None):
 
     Returns
     -------
-    opening : uint8 array
+    out : array, same shape and type as `image`
         The result of the morphological white top hat.
 
     Examples
@@ -281,16 +395,18 @@ def white_tophat(image, selem=None, out=None):
            [0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-
-    if image is out:
-        raise NotImplementedError("Cannot perform white top hat in place.")
-
-    out = opening(image, selem, out=out)
-    out = image - out
+    selem = np.array(selem)
+    if out is image:
+        opened = opening(image, selem)
+        out -= opened
+        return out
+    elif out is None:
+        out = np.empty_like(image)
+    out = nd.white_tophat(image, footprint=selem, output=out)
     return out
 
 
-@default_fallback
+@default_selem
 def black_tophat(image, selem=None, out=None):
     """Return black top hat of an image.
 
@@ -312,7 +428,7 @@ def black_tophat(image, selem=None, out=None):
 
     Returns
     -------
-    opening : uint8 array
+    opening : array, same shape and type as `image`
        The result of the black top filter.
 
     Examples
@@ -333,10 +449,10 @@ def black_tophat(image, selem=None, out=None):
            [0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-
-    if image is out:
-        raise NotImplementedError("Cannot perform white top hat in place.")
-
+    if out is image:
+        original = image.copy()
+    else:
+        original = image
     out = closing(image, selem, out=out)
-    out = out - image
+    out -= original
     return out
diff --git a/skimage/morphology/misc.py b/skimage/morphology/misc.py
index 55b279a2..663e1e02 100644
--- a/skimage/morphology/misc.py
+++ b/skimage/morphology/misc.py
@@ -15,11 +15,8 @@ funcs = ('binary_erosion', 'binary_dilation', 'binary_opening',
 skimage2ndimage.update(dict((x, x) for x in funcs))
 
 
-def default_fallback(func):
-    """Decorator to fall back on ndimage for images with more than 2 dimensions
-
-    Decorator also provides a default structuring element, `selem`, with the
-    appropriate dimensionality if none is specified.
+def default_selem(func):
+    """Decorator to add a default structuring element to morphology functions.
 
     Parameters
     ----------
@@ -30,25 +27,14 @@ def default_fallback(func):
     Returns
     -------
     func_out : function
-        If the image dimensionality is greater than 2D, the ndimage
-        function is returned, otherwise skimage function is used.
+        The function, using a default structuring element of same dimension
+        as the input image with connectivity 1.
     """
     @functools.wraps(func)
-    def func_out(image, selem=None, out=None, **kwargs):
-        # Default structure element
+    def func_out(image, selem=None, *args, **kwargs):
         if selem is None:
             selem = _default_selem(image.ndim)
-
-        # If image has more than 2 dimensions, use scipy.ndimage
-        if image.ndim > 2:
-            function = getattr(nd, skimage2ndimage[func.__name__])
-            try:
-                return function(image, footprint=selem, output=out, **kwargs)
-            except TypeError:
-                # nd.binary_* take structure instead of footprint
-                return function(image, structure=selem, output=out, **kwargs)
-        else:
-            return func(image, selem=selem, out=out, **kwargs)
+        return func(image, selem=selem, *args, **kwargs)
 
     return func_out
 
diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py
index de1a8794..dbbcad8b 100644
--- a/skimage/morphology/setup.py
+++ b/skimage/morphology/setup.py
@@ -12,14 +12,11 @@ def configuration(parent_package='', top_path=None):
     config = Configuration('morphology', parent_package, top_path)
     config.add_data_dir('tests')
 
-    cython(['cmorph.pyx'], working_path=base_path)
     cython(['_watershed.pyx'], working_path=base_path)
     cython(['_skeletonize_cy.pyx'], working_path=base_path)
     cython(['_convex_hull.pyx'], working_path=base_path)
     cython(['_greyreconstruct.pyx'], working_path=base_path)
 
-    config.add_extension('cmorph', sources=['cmorph.c'],
-                         include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_watershed', sources=['_watershed.c'],
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_skeletonize_cy', sources=['_skeletonize_cy.c'],
diff --git a/skimage/morphology/tests/test_binary.py b/skimage/morphology/tests/test_binary.py
index d52f92bb..6d8097b4 100644
--- a/skimage/morphology/tests/test_binary.py
+++ b/skimage/morphology/tests/test_binary.py
@@ -4,7 +4,6 @@ from numpy import testing
 from skimage import data, color
 from skimage.util import img_as_bool
 from skimage.morphology import binary, grey, selem
-from skimage._shared._warnings import expected_warnings
 from scipy import ndimage
 
 
@@ -15,50 +14,44 @@ bw_img = img > 100
 def test_non_square_image():
     strel = selem.square(3)
     binary_res = binary.binary_erosion(bw_img[:100, :200], strel)
-    with expected_warnings(['precision loss']):
-        grey_res = img_as_bool(grey.erosion(bw_img[:100, :200], strel))
+    grey_res = img_as_bool(grey.erosion(bw_img[:100, :200], strel))
     testing.assert_array_equal(binary_res, grey_res)
 
 
 def test_binary_erosion():
     strel = selem.square(3)
     binary_res = binary.binary_erosion(bw_img, strel)
-    with expected_warnings(['precision loss']):
-        grey_res = img_as_bool(grey.erosion(bw_img, strel))
+    grey_res = img_as_bool(grey.erosion(bw_img, strel))
     testing.assert_array_equal(binary_res, grey_res)
 
 
 def test_binary_dilation():
     strel = selem.square(3)
     binary_res = binary.binary_dilation(bw_img, strel)
-    with expected_warnings(['precision loss']):
-        grey_res = img_as_bool(grey.dilation(bw_img, strel))
+    grey_res = img_as_bool(grey.dilation(bw_img, strel))
     testing.assert_array_equal(binary_res, grey_res)
 
 
 def test_binary_closing():
     strel = selem.square(3)
     binary_res = binary.binary_closing(bw_img, strel)
-    with expected_warnings(['precision loss']):
-        grey_res = img_as_bool(grey.closing(bw_img, strel))
+    grey_res = img_as_bool(grey.closing(bw_img, strel))
     testing.assert_array_equal(binary_res, grey_res)
 
 
 def test_binary_opening():
     strel = selem.square(3)
     binary_res = binary.binary_opening(bw_img, strel)
-    with expected_warnings(['precision loss']):
-        grey_res = img_as_bool(grey.opening(bw_img, strel))
+    grey_res = img_as_bool(grey.opening(bw_img, strel))
     testing.assert_array_equal(binary_res, grey_res)
 
 
 def test_selem_overflow():
     strel = np.ones((17, 17), dtype=np.uint8)
-    img = np.zeros((20, 20))
-    img[2:19, 2:19] = 1
+    img = np.zeros((20, 20), dtype=bool)
+    img[2:19, 2:19] = True
     binary_res = binary.binary_erosion(img, strel)
-    with expected_warnings(['precision loss']):
-        grey_res = img_as_bool(grey.erosion(img, strel))
+    grey_res = img_as_bool(grey.erosion(img, strel))
     testing.assert_array_equal(binary_res, grey_res)
 
 
diff --git a/skimage/morphology/tests/test_grey.py b/skimage/morphology/tests/test_grey.py
index 911e0c71..576404f2 100644
--- a/skimage/morphology/tests/test_grey.py
+++ b/skimage/morphology/tests/test_grey.py
@@ -5,7 +5,7 @@ from numpy import testing
 from scipy import ndimage
 
 import skimage
-from skimage import data_dir
+from skimage import data_dir, img_as_uint
 from skimage.morphology import grey, selem
 from skimage._shared._warnings import expected_warnings
 
@@ -208,51 +208,51 @@ def test_2d_ndimage_equivalence():
     testing.assert_array_equal(opened, ndimage_opened)
     testing.assert_array_equal(closed, ndimage_closed)
 
-class TestDTypes():
+# float test images
+im = np.array([[ 0.55,  0.72,  0.6 ,  0.54,  0.42],
+               [ 0.65,  0.44,  0.89,  0.96,  0.38],
+               [ 0.79,  0.53,  0.57,  0.93,  0.07],
+               [ 0.09,  0.02,  0.83,  0.78,  0.87],
+               [ 0.98,  0.8 ,  0.46,  0.78,  0.12]])
 
-    def setUp(self):
-        k = 5
-        arrname = '%03i' % k
+eroded = np.array([[ 0.55,  0.44,  0.54,  0.42,  0.38],
+                   [ 0.44,  0.44,  0.44,  0.38,  0.07],
+                   [ 0.09,  0.02,  0.53,  0.07,  0.07],
+                   [ 0.02,  0.02,  0.02,  0.78,  0.07],
+                   [ 0.09,  0.02,  0.46,  0.12,  0.12]])
 
-        self.disk = selem.disk(k)
+dilated = np.array([[ 0.72,  0.72,  0.89,  0.96,  0.54],
+                    [ 0.79,  0.89,  0.96,  0.96,  0.96],
+                    [ 0.79,  0.79,  0.93,  0.96,  0.93],
+                    [ 0.98,  0.83,  0.83,  0.93,  0.87],
+                    [ 0.98,  0.98,  0.83,  0.78,  0.87]])
 
-        fname_opening = os.path.join(data_dir, "disk-open-matlab-output.npz")
-        self.expected_opening = np.load(fname_opening)[arrname]
+opened = np.array([[ 0.55,  0.55,  0.54,  0.54,  0.42],
+                   [ 0.55,  0.44,  0.54,  0.44,  0.38],
+                   [ 0.44,  0.53,  0.53,  0.78,  0.07],
+                   [ 0.09,  0.02,  0.78,  0.78,  0.78],
+                   [ 0.09,  0.46,  0.46,  0.78,  0.12]])
 
-        fname_closing = os.path.join(data_dir, "disk-close-matlab-output.npz")
-        self.expected_closing = np.load(fname_closing)[arrname]
+closed = np.array([[ 0.72,  0.72,  0.72,  0.54,  0.54],
+                   [ 0.72,  0.72,  0.89,  0.96,  0.54],
+                   [ 0.79,  0.79,  0.79,  0.93,  0.87],
+                   [ 0.79,  0.79,  0.83,  0.78,  0.87],
+                   [ 0.98,  0.83,  0.78,  0.78,  0.78]])
 
-    def _test_image(self, image):
-        with expected_warnings(['precision loss']):
-            result_opening = grey.opening(image, self.disk)
-        testing.assert_equal(result_opening, self.expected_opening)
-
-        with expected_warnings(['precision loss']):
-            result_closing = grey.closing(image, self.disk)
-        testing.assert_equal(result_closing, self.expected_closing)
-
-    def test_float(self):
-        image = skimage.img_as_float(lena)
-        self._test_image(image)
-
-    @testing.decorators.skipif(True)
-    def test_int(self):
-        image = skimage.img_as_int(lena)
-        self._test_image(image)
-
-    def test_uint(self):
-        image = skimage.img_as_uint(lena)
-        self._test_image(image)
+def test_float():
+    np.testing.assert_allclose(grey.erosion(im), eroded)
+    np.testing.assert_allclose(grey.dilation(im), dilated)
+    np.testing.assert_allclose(grey.opening(im), opened)
+    np.testing.assert_allclose(grey.closing(im), closed)
 
 
-def test_inplace():
-    selem = np.ones((3, 3))
-    image = np.zeros((5, 5))
-    out = image
-
-    for f in (grey.erosion, grey.dilation,
-              grey.white_tophat, grey.black_tophat):
-        testing.assert_raises(NotImplementedError, f, image, selem, out=out)
+def test_uint16():
+    im16, eroded16, dilated16, opened16, closed16 = (
+                    map(img_as_uint, [im, eroded, dilated, opened, closed]))
+    np.testing.assert_allclose(grey.erosion(im16), eroded16)
+    np.testing.assert_allclose(grey.dilation(im16), dilated16)
+    np.testing.assert_allclose(grey.opening(im16), opened16)
+    np.testing.assert_allclose(grey.closing(im16), closed16)
 
 
 def test_discontiguous_out_array():
diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py
index 8e4ef32e..1975d879 100644
--- a/skimage/segmentation/boundaries.py
+++ b/skimage/segmentation/boundaries.py
@@ -1,41 +1,218 @@
+from __future__ import division
+
 import numpy as np
-from ..morphology import dilation, square
-from ..util import img_as_float
+from scipy import ndimage as nd
+from ..morphology import dilation, erosion, square
+from ..util import img_as_float, view_as_windows, pad
 from ..color import gray2rgb
-from .._shared.utils import deprecated
 
 
-def find_boundaries(label_img):
-    """Return bool array where boundaries between labeled regions are True."""
-    boundaries = np.zeros(label_img.shape, dtype=np.bool)
-    boundaries[1:, :] += label_img[1:, :] != label_img[:-1, :]
-    boundaries[:, 1:] += label_img[:, 1:] != label_img[:, :-1]
+def _find_boundaries_subpixel(label_img):
+    """See ``find_boundaries(..., mode='subpixel')``.
+
+    Notes
+    -----
+    This function puts in an empty row and column between each *actual*
+    row and column of the image, for a corresponding shape of $2s - 1$
+    for every image dimension of size $s$. These "interstitial" rows
+    and columns are filled as ``True`` if they separate two labels in
+    `label_img`, ``False`` otherwise.
+
+    I used ``view_as_windows`` to get the neighborhood of each pixel.
+    Then I check whether there are two labels or more in that
+    neighborhood.
+    """
+    ndim = label_img.ndim
+    max_label = np.iinfo(label_img.dtype).max
+
+    label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape],
+                                  label_img.dtype)
+    pixels = [slice(None, None, 2)] * ndim
+    label_img_expanded[pixels] = label_img
+
+    edges = np.ones(label_img_expanded.shape, dtype=bool)
+    edges[pixels] = False
+    label_img_expanded[edges] = max_label
+    windows = view_as_windows(pad(label_img_expanded, 1,
+                                  mode='constant', constant_values=0),
+                              (3,) * ndim)
+
+    boundaries = np.zeros_like(edges)
+    for index in np.ndindex(label_img_expanded.shape):
+        if edges[index]:
+            values = np.unique(windows[index].ravel())
+            if len(values) > 2: # single value and max_label
+                boundaries[index] = True
     return boundaries
 
 
+def find_boundaries(label_img, connectivity=1, mode='thick', background=0):
+    """Return bool array where boundaries between labeled regions are True.
+
+    Parameters
+    ----------
+    label_img : array of int
+        An array in which different regions are labeled with different
+        integers.
+    connectivity: int in {1, ..., `label_img.ndim`}, optional
+        A pixel is considered a boundary pixel if any of its neighbors
+        has a different label. `connectivity` controls which pixels are
+        considered neighbors. A connectivity of 1 (default) means
+        pixels sharing an edge (in 2D) or a face (in 3D) will be
+        considered neighbors. A connectivity of `label_img.ndim` means
+        pixels sharing a corner will be considered neighbors.
+    mode: string in {'thick', 'inner', 'outer', 'subpixel'}
+        How to mark the boundaries:
+
+        - thick: any pixel not completely surrounded by pixels of the
+        same label (defined by `connectivity`) is marked as a boundary.
+        This results in boundaries that are 2 pixels thick.
+        - inner: outline the pixels *just inside* of objects, leaving
+        background pixels untouched.
+        - outer: outline pixels in the background around object
+        boundaries. When two objects touch, their boundary is also
+        marked.
+        - subpixel: return a doubled image, with pixels *between* the
+        original pixels marked as boundary where appropriate.
+    background: int, optional
+        For modes 'inner' and 'outer', a definition of a background
+        label is required. See `mode` for descriptions of these two.
+
+    Returns
+    -------
+    boundaries : array of bool, same shape as `label_img`
+        A bool image where ``True`` represents a boundary pixel. For
+        `mode` equal to 'subpixel', ``boundaries.shape[i]`` is equal
+        to ``2 * label_img.shape[i] - 1`` for all ``i`` (a pixel is
+        inserted in between all other pairs of pixels).
+
+    Examples
+    --------
+    >>> labels = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ...                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ...                    [0, 0, 0, 0, 0, 5, 5, 5, 0, 0],
+    ...                    [0, 0, 1, 1, 1, 5, 5, 5, 0, 0],
+    ...                    [0, 0, 1, 1, 1, 5, 5, 5, 0, 0],
+    ...                    [0, 0, 1, 1, 1, 5, 5, 5, 0, 0],
+    ...                    [0, 0, 0, 0, 0, 5, 5, 5, 0, 0],
+    ...                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    ...                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8)
+    >>> find_boundaries(labels, mode='thick').astype(np.uint8)
+    array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 1, 1, 0],
+           [0, 1, 1, 1, 1, 1, 0, 1, 1, 0],
+           [0, 1, 1, 0, 1, 1, 0, 1, 1, 0],
+           [0, 1, 1, 1, 1, 1, 0, 1, 1, 0],
+           [0, 0, 1, 1, 1, 1, 1, 1, 1, 0],
+           [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+    >>> find_boundaries(labels, mode='inner').astype(np.uint8)
+    array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 0, 1, 0, 0],
+           [0, 0, 1, 0, 1, 1, 0, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 0, 1, 0, 0],
+           [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+    >>> find_boundaries(labels, mode='outer').astype(np.uint8)
+    array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 0, 0, 1, 0],
+           [0, 1, 0, 0, 1, 1, 0, 0, 1, 0],
+           [0, 1, 0, 0, 1, 1, 0, 0, 1, 0],
+           [0, 1, 0, 0, 1, 1, 0, 0, 1, 0],
+           [0, 0, 1, 1, 1, 1, 0, 0, 1, 0],
+           [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+    >>> labels_small = labels[::2, ::3]
+    >>> labels_small
+    array([[0, 0, 0, 0],
+           [0, 0, 5, 0],
+           [0, 1, 5, 0],
+           [0, 0, 5, 0],
+           [0, 0, 0, 0]], dtype=uint8)
+    >>> find_boundaries(labels_small, mode='subpixel').astype(np.uint8)
+    array([[0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 1, 1, 1, 0],
+           [0, 0, 0, 1, 0, 1, 0],
+           [0, 1, 1, 1, 0, 1, 0],
+           [0, 1, 0, 1, 0, 1, 0],
+           [0, 1, 1, 1, 0, 1, 0],
+           [0, 0, 0, 1, 0, 1, 0],
+           [0, 0, 0, 1, 1, 1, 0],
+           [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+    """
+    ndim = label_img.ndim
+    selem = nd.generate_binary_structure(ndim, connectivity)
+    if mode != 'subpixel':
+        boundaries = dilation(label_img, selem) != erosion(label_img, selem)
+        if mode == 'inner':
+            foreground_image = (label_img != background)
+            boundaries &= foreground_image
+        elif mode == 'outer':
+            max_label = np.iinfo(label_img.dtype).max
+            background_image = (label_img == background)
+            selem = nd.generate_binary_structure(ndim, ndim)
+            inverted_background = np.array(label_img, copy=True)
+            inverted_background[background_image] = max_label
+            adjacent_objects = ((dilation(label_img, selem) !=
+                                 erosion(inverted_background, selem)) &
+                                ~background_image)
+            boundaries &= (background_image | adjacent_objects)
+        return boundaries
+    else:
+        boundaries = _find_boundaries_subpixel(label_img)
+        return boundaries
+
+
 def mark_boundaries(image, label_img, color=(1, 1, 0),
-                    outline_color=(0, 0, 0)):
+                    outline_color=None, mode='outer', background_label=0):
     """Return image with boundaries between labeled regions highlighted.
 
     Parameters
     ----------
     image : (M, N[, 3]) array
         Grayscale or RGB image.
-    label_img : (M, N) array
+    label_img : (M, N) array of int
         Label array where regions are marked by different integer values.
-    color : length-3 sequence
+    color : length-3 sequence, optional
         RGB color of boundaries in the output image.
-    outline_color : length-3 sequence
+    outline_color : length-3 sequence, optional
         RGB color surrounding boundaries in the output image. If None, no
         outline is drawn.
-    """
-    if image.ndim == 2:
-        image = gray2rgb(image)
-    image = img_as_float(image, force_copy=True)
+    mode : string in {'thick', 'inner', 'outer', 'subpixel'}, optional
+        The mode for finding boundaries.
+    background_label : int, optional
+        Which label to consider background (this is only useful for
+        modes ``inner`` and ``outer``).
 
-    boundaries = find_boundaries(label_img)
+    Returns
+    -------
+    marked : (M, N, 3) array of float
+        An image in which the boundaries between labels are
+        superimposed on the original image.
+
+    See Also
+    --------
+    ``find_boundaries``.
+    """
+    marked = img_as_float(image, force_copy=True)
+    if marked.ndim == 2:
+        marked = gray2rgb(marked)
+    if mode == 'subpixel':
+        # Here, we want to interpose an extra line of pixels between
+        # each original line - except for the last axis which holds
+        # the RGB information. ``nd.zoom`` then performs the (cubic)
+        # interpolation, filling in the values of the interposed pixels
+        marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1],
+                         mode='reflect')
+    boundaries = find_boundaries(label_img, mode=mode,
+                                 background=background_label)
     if outline_color is not None:
-        outer_boundaries = dilation(boundaries.astype(np.uint8), square(2))
-        image[outer_boundaries != 0, :] = np.array(outline_color)
-    image[boundaries, :] = np.array(color)
-    return image
+        outlines = dilation(boundaries, square(3))
+        marked[outlines] = outline_color
+    marked[boundaries] = color
+    return marked
diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py
index 2fff52f8..e6e401ac 100644
--- a/skimage/segmentation/tests/test_boundaries.py
+++ b/skimage/segmentation/tests/test_boundaries.py
@@ -1,19 +1,22 @@
 import numpy as np
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_equal, assert_allclose
 from skimage.segmentation import find_boundaries, mark_boundaries
 
 
+white = (1, 1, 1)
+
+
 def test_find_boundaries():
-    image = np.zeros((10, 10))
+    image = np.zeros((10, 10), dtype=np.uint8)
     image[2:7, 2:7] = 1
 
     ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                    [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
+                    [0, 0, 1, 1, 1, 1, 1, 0, 0, 0],
+                    [0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
                     [0, 0, 1, 1, 1, 1, 1, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
@@ -24,36 +27,63 @@ def test_find_boundaries():
 
 def test_mark_boundaries():
     image = np.zeros((10, 10))
-    label_image = np.zeros((10, 10))
+    label_image = np.zeros((10, 10), dtype=np.uint8)
     label_image[2:7, 2:7] = 1
 
     ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                    [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
-                    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
+                    [0, 0, 1, 1, 1, 1, 1, 0, 0, 0],
+                    [0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 1, 1, 0, 0],
+                    [0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
                     [0, 0, 1, 1, 1, 1, 1, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
-    result = mark_boundaries(image, label_image, color=(1, 1, 1)).mean(axis=2)
+
+    marked = mark_boundaries(image, label_image, color=white, mode='thick')
+    result = np.mean(marked, axis=-1)
     assert_array_equal(result, ref)
 
-    ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                    [0, 0, 1, 1, 1, 1, 1, 1, 2, 0],
-                    [0, 0, 1, 2, 2, 2, 2, 1, 2, 0],
-                    [0, 0, 1, 2, 0, 0, 0, 1, 2, 0],
-                    [0, 0, 1, 2, 0, 0, 0, 1, 2, 0],
-                    [0, 0, 1, 2, 0, 0, 0, 1, 2, 0],
-                    [0, 0, 1, 1, 1, 1, 1, 2, 2, 0],
-                    [0, 0, 2, 2, 2, 2, 2, 2, 0, 0],
+    ref = np.array([[0, 2, 2, 2, 2, 2, 2, 2, 0, 0],
+                    [2, 2, 1, 1, 1, 1, 1, 2, 2, 0],
+                    [2, 1, 1, 1, 1, 1, 1, 1, 2, 0],
+                    [2, 1, 1, 2, 2, 2, 1, 1, 2, 0],
+                    [2, 1, 1, 2, 0, 2, 1, 1, 2, 0],
+                    [2, 1, 1, 2, 2, 2, 1, 1, 2, 0],
+                    [2, 1, 1, 1, 1, 1, 1, 1, 2, 0],
+                    [2, 2, 1, 1, 1, 1, 1, 2, 2, 0],
+                    [0, 2, 2, 2, 2, 2, 2, 2, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
-    result = mark_boundaries(image, label_image, color=(1, 1, 1),
-                             outline_color=(2, 2, 2)).mean(axis=2)
+    marked = mark_boundaries(image, label_image, color=white,
+                             outline_color=(2, 2, 2), mode='thick')
+    result = np.mean(marked, axis=-1)
     assert_array_equal(result, ref)
 
 
+def test_mark_boundaries_subpixel():
+    labels = np.array([[0, 0, 0, 0],
+                       [0, 0, 5, 0],
+                       [0, 1, 5, 0],
+                       [0, 0, 5, 0],
+                       [0, 0, 0, 0]], dtype=np.uint8)
+    np.random.seed(0)
+    image = np.round(np.random.rand(*labels.shape), 2)
+    marked = mark_boundaries(image, labels, color=white, mode='subpixel')
+    marked_proj = np.round(np.mean(marked, axis=-1), 2)
+
+    ref_result = np.array(
+            [[ 0.55,  0.63,  0.72,  0.69,  0.6 ,  0.55,  0.54],
+             [ 0.45,  0.58,  0.72,  1.  ,  1.  ,  1.  ,  0.69],
+             [ 0.42,  0.54,  0.65,  1.  ,  0.44,  1.  ,  0.89],
+             [ 0.69,  1.  ,  1.  ,  1.  ,  0.69,  1.  ,  0.83],
+             [ 0.96,  1.  ,  0.38,  1.  ,  0.79,  1.  ,  0.53],
+             [ 0.89,  1.  ,  1.  ,  1.  ,  0.38,  1.  ,  0.16],
+             [ 0.57,  0.78,  0.93,  1.  ,  0.07,  1.  ,  0.09],
+             [ 0.2 ,  0.52,  0.92,  1.  ,  1.  ,  1.  ,  0.54],
+             [ 0.02,  0.35,  0.83,  0.9 ,  0.78,  0.81,  0.87]])
+    assert_allclose(marked_proj, ref_result, atol=0.01)
+
+
 if __name__ == "__main__":
     np.testing.run_module_suite()
diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index 5577e46b..4739bd5b 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -3,7 +3,7 @@ from .dtype import (img_as_float, img_as_int, img_as_uint, img_as_ubyte,
 from .shape import view_as_blocks, view_as_windows
 from .noise import random_noise
 
-from .arraypad import pad
+from .arraypad import pad, crop
 from ._regular_grid import regular_grid
 from .unique import unique_rows
 
@@ -17,6 +17,7 @@ __all__ = ['img_as_float',
            'view_as_blocks',
            'view_as_windows',
            'pad',
+           'crop',
            'random_noise',
            'regular_grid',
            'unique_rows']
diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py
index d616102d..f8104df5 100644
--- a/skimage/util/arraypad.py
+++ b/skimage/util/arraypad.py
@@ -8,7 +8,7 @@ from __future__ import division, absolute_import, print_function
 import numpy as np
 
 
-__all__ = ['pad']
+__all__ = ['pad', 'crop']
 
 
 ###############################################################################
@@ -1496,3 +1496,43 @@ def pad(array, pad_width, mode=None, **kwargs):
             newmat = _pad_wrap(newmat, (pad_before, pad_after), axis)
 
     return newmat
+
+
+def crop(ar, crop_width, copy=False, order='K'):
+    """Crop array `ar` by `crop_width` along each dimension.
+
+    Parameters
+    ----------
+    ar : array-like of rank N
+        Input array.
+    crop_width : {sequence, int}
+        Number of values to remove from the edges of each axis.
+        ``((before_1, after_1),`` ... ``(before_N, after_N))`` specifies
+        unique crop widths at the start and end of each axis.
+        ``((before, after),)`` specifies a fixed start and end crop
+        for every axis.
+        ``(n,)`` or ``n`` for integer ``n`` is a shortcut for
+        before = after = ``n`` for all axes.
+    copy : bool, optional
+        Ensure the returned array is a contiguous copy. Normally, a crop
+        operation will return a discontiguous view of the underlying
+        input array. Passing ``copy=True`` will result in a contiguous
+        copy.
+    order : {'C', 'F', 'A', 'K'}, optional
+        If ``copy==True``, control the memory layout of the copy. See
+        ``np.copy``.
+
+    Returns
+    -------
+    cropped : array
+        The cropped array. If ``copy=False`` (default), this is a sliced
+        view of the input array.
+    """
+    ar = np.array(ar, copy=False)
+    crops = _validate_lengths(ar, crop_width)
+    slices = [slice(a, ar.shape[i] - b) for i, (a, b) in enumerate(crops)]
+    if copy:
+        cropped = np.array(ar[slices], order=order, copy=True)
+    else:
+        cropped = ar[slices]
+    return cropped
diff --git a/skimage/util/tests/test_arraypad.py b/skimage/util/tests/test_arraypad.py
index eb8554b1..d6f993f9 100644
--- a/skimage/util/tests/test_arraypad.py
+++ b/skimage/util/tests/test_arraypad.py
@@ -5,8 +5,8 @@ from __future__ import division, absolute_import, print_function
 
 import numpy as np
 from numpy.testing import (assert_array_equal, assert_raises, assert_allclose,
-                           TestCase)
-from skimage.util import pad
+                           assert_equal, TestCase)
+from skimage.util import pad, crop
 
 
 class TestConditionalShortcuts(TestCase):
@@ -1043,5 +1043,49 @@ class TypeError1(TestCase):
                       **kwargs)
 
 
+def test_multi_crop():
+    arr = np.arange(45).reshape(9, 5)
+    out = crop(arr, ((1, 2), (2, 1)))
+    assert_array_equal(out[0], [7, 8])
+    assert_array_equal(out[-1], [32, 33])
+    assert_equal(out.shape, (6, 2))
+
+
+def test_pair_crop():
+    arr = np.arange(45).reshape(9, 5)
+    out = crop(arr, (1, 2))
+    assert_array_equal(out[0], [6, 7])
+    assert_array_equal(out[-1], [31, 32])
+    assert_equal(out.shape, (6, 2))
+
+
+def test_int_crop():
+    arr = np.arange(45).reshape(9, 5)
+    out = crop(arr, 1)
+    assert_array_equal(out[0], [6, 7, 8])
+    assert_array_equal(out[-1], [36, 37, 38])
+    assert_equal(out.shape, (7, 3))
+
+
+def test_copy_crop():
+    arr = np.arange(45).reshape(9, 5)
+    out0 = crop(arr, 1, copy=True)
+    assert out0.flags.c_contiguous
+    out0[0, 0] = 100
+    assert not np.any(arr == 100)
+    assert not np.may_share_memory(arr, out0)
+
+    out1 = crop(arr, 1)
+    out1[0, 0] = 100
+    assert arr[1, 1] == 100
+    assert np.may_share_memory(arr, out1)
+
+
+def test_zero_crop():
+    arr = np.arange(45).reshape(9, 5)
+    out = crop(arr, 0)
+    assert out.shape == (9, 5)
+
+
 if __name__ == "__main__":
     np.testing.run_module_suite()