diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index d57b6015009..92ac35b3128 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -202,4 +202,7 @@
   Extended the image labelling implementation so it also works on 3D images.
 
 - Salvatore Scaramuzzino
-  RectTool example
\ No newline at end of file
+  RectTool example
+
+- Kevin Keraudren
+  Fix and test for feature.peak_local_max
diff --git a/bento.info b/bento.info
index 6a8a0000d17..a36a4ebb426 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/__init__.py b/skimage/__init__.py
index fb38f87a7f2..303ea34fe5d 100644
--- a/skimage/__init__.py
+++ b/skimage/__init__.py
@@ -60,7 +60,6 @@
 import imp as _imp
 import functools as _functools
 import warnings as _warnings
-from skimage._shared.utils import deprecated as _deprecated
 
 pkg_dir = _osp.abspath(_osp.dirname(__file__))
 data_dir = _osp.join(pkg_dir, 'data')
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index 615d9937669..649f69f27a6 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -19,8 +19,8 @@ def _coords_inside_image(rr, cc, shape, val=None):
     rr, cc : (N,) ndarray of int
         Indices of pixels.
     shape : tuple
-        Image shape which is used to determine maximum extents of output pixel
-        coordinates.
+        Image shape which is used to determine the maximum extent of output
+        pixel coordinates.
     val : ndarray of float, optional
         Values of pixels at coordinates [rr, cc].
 
@@ -223,9 +223,9 @@ def polygon(y, x, shape=None):
     x : (N,) ndarray
         X-coordinates of vertices of polygon.
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
-        coordinates. This is useful for polygons which exceed the image size.
-        By default the full extents of the polygon are used.
+        Image shape which is used to determine the maximum extent of output
+        pixel coordinates. This is useful for polygons which exceed the image
+        size. By default the full extent of the polygon are used.
 
     Returns
     -------
@@ -303,9 +303,9 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
         bresenham : Bresenham method (default)
         andres : Andres method
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
+        Image shape which is used to determine the maximum extent of output pixel
         coordinates. This is useful for circles which exceed the image size.
-        By default the full extents of the polygon are used.
+        By default the full extent of the circle are used.
 
     Returns
     -------
@@ -411,9 +411,9 @@ def circle_perimeter_aa(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
     radius: int
         Radius of circle.
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
+        Image shape which is used to determine the maximum extent of output pixel
         coordinates. This is useful for circles which exceed the image size.
-        By default the full extents of the polygon are used.
+        By default the full extent of the circle are used.
 
     Returns
     -------
@@ -499,9 +499,9 @@ def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius,
     orientation : double, optional (default 0)
         Major axis orientation in clockwise direction as radians.
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
+        Image shape which is used to determine the maximum extent of output pixel
         coordinates. This is useful for ellipses which exceed the image size.
-        By default the full extents of the polygon are used.
+        By default the full extent of the ellipse are used.
 
     Returns
     -------
@@ -774,9 +774,9 @@ def bezier_curve(Py_ssize_t y0, Py_ssize_t x0,
     weight : double
         Middle control point weight, it describes the line tension.
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
-        coordinates. This is useful for curves which exceed the image size.
-        By default the full extents of the polygon are used.
+        Image shape which is used to determine the maximum extent of output
+        pixel coordinates. This is useful for curves which exceed the image
+        size. By default the full extent of the curve are used.
 
     Returns
     -------
diff --git a/skimage/draw/draw.py b/skimage/draw/draw.py
index b6bc5588baf..be59e268afd 100644
--- a/skimage/draw/draw.py
+++ b/skimage/draw/draw.py
@@ -22,9 +22,9 @@ def ellipse(cy, cx, yradius, xradius, shape=None):
     yradius, xradius : double
         Minor and major semi-axes. ``(x/xradius)**2 + (y/yradius)**2 = 1``.
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
+        Image shape which is used to determine the maximum extent of output pixel
         coordinates. This is useful for ellipses which exceed the image size.
-        By default the full extents of the ellipse are used.
+        By default the full extent of the ellipse are used.
 
     Returns
     -------
@@ -85,9 +85,9 @@ def circle(cy, cx, radius, shape=None):
     radius: double
         Radius of circle.
     shape : tuple, optional
-        Image shape which is used to determine maximum extents of output pixel
+        Image shape which is used to determine the maximum extent of output pixel
         coordinates. This is useful for circles which exceed the image size.
-        By default the full extents of the circle are used.
+        By default the full extent of the circle are used.
 
     Returns
     -------
diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index e70f6df16f4..599a6e581cd 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -146,10 +146,10 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     peak_threshold = max(np.max(image.ravel()) * threshold_rel, threshold_abs)
 
     # get coordinates of peaks
-    coordinates = np.transpose((image > peak_threshold).nonzero())
+    coordinates = np.argwhere(image > peak_threshold)
 
     if coordinates.shape[0] > num_peaks:
-        intensities = image[coordinates[:, 0], coordinates[:, 1]]
+        intensities = image.flat[np.ravel_multi_index(coordinates.transpose(),image.shape)]
         idx_maxsort = np.argsort(intensities)[::-1]
         coordinates = coordinates[idx_maxsort][:num_peaks]
 
diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py
index ec130c663ff..09e524953df 100644
--- a/skimage/feature/tests/test_peak.py
+++ b/skimage/feature/tests/test_peak.py
@@ -82,6 +82,15 @@ def test_num_peaks():
     assert (3, 5) in peaks_limited
 
 
+def test_num_peaks3D():
+    # Issue 1354: the old code only hold for 2D arrays
+    # and this code would die with IndexError
+    image = np.zeros((10, 10, 100))
+    image[5,5,::5] = np.arange(20)
+    peaks_limited = peak.peak_local_max(image, min_distance=1, num_peaks=2)
+    assert len(peaks_limited) == 2
+    
+
 def test_reorder_labels():
     image = np.random.uniform(size=(40, 60))
     i, j = np.mgrid[0:40, 0:60]
diff --git a/skimage/filters/rank/README.rst b/skimage/filters/rank/README.rst
index e5c5a9ad26e..6e063ec2737 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 64c3ca979e0..4862843253d 100644
--- a/skimage/filters/rank/tests/test_rank.py
+++ b/skimage/filters/rank/tests/test_rank.py
@@ -5,7 +5,7 @@
 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/measure/_regionprops.py b/skimage/measure/_regionprops.py
index 66a8734f055..20e7e590b14 100644
--- a/skimage/measure/_regionprops.py
+++ b/skimage/measure/_regionprops.py
@@ -115,21 +115,21 @@ def __init__(self, slice, label, label_image, intensity_image,
         self._intensity_image = intensity_image
         self._cache_active = cache_active
 
-    @_cached_property
+    @property
     def area(self):
         return self.moments[0, 0]
 
-    @_cached_property
+    @property
     def bbox(self):
         return (self._slice[0].start, self._slice[1].start,
                 self._slice[0].stop, self._slice[1].stop)
 
-    @_cached_property
+    @property
     def centroid(self):
         row, col = self.local_centroid
         return row + self._slice[0].start, col + self._slice[1].start
 
-    @_cached_property
+    @property
     def convex_area(self):
         return np.sum(self.convex_image)
 
@@ -138,35 +138,35 @@ def convex_image(self):
         from ..morphology.convex_hull import convex_hull_image
         return convex_hull_image(self.image)
 
-    @_cached_property
+    @property
     def coords(self):
         rr, cc = np.nonzero(self.image)
         return np.vstack((rr + self._slice[0].start,
                           cc + self._slice[1].start)).T
 
-    @_cached_property
+    @property
     def eccentricity(self):
         l1, l2 = self.inertia_tensor_eigvals
         if l1 == 0:
             return 0
         return sqrt(1 - l2 / l1)
 
-    @_cached_property
+    @property
     def equivalent_diameter(self):
         return sqrt(4 * self.moments[0, 0] / PI)
 
-    @_cached_property
+    @property
     def euler_number(self):
         euler_array = self.filled_image != self.image
         _, num = label(euler_array, neighbors=8, return_num=True)
         return -num + 1
 
-    @_cached_property
+    @property
     def extent(self):
         rows, cols = self.image.shape
         return self.moments[0, 0] / (rows * cols)
 
-    @_cached_property
+    @property
     def filled_area(self):
         return np.sum(self.filled_image)
 
@@ -200,35 +200,35 @@ def intensity_image(self):
             raise AttributeError('No intensity image specified.')
         return self._intensity_image[self._slice] * self.image
 
-    @_cached_property
+    @property
     def _intensity_image_double(self):
         return self.intensity_image.astype(np.double)
 
-    @_cached_property
+    @property
     def local_centroid(self):
         m = self.moments
         row = m[0, 1] / m[0, 0]
         col = m[1, 0] / m[0, 0]
         return row, col
 
-    @_cached_property
+    @property
     def max_intensity(self):
         return np.max(self.intensity_image[self.image])
 
-    @_cached_property
+    @property
     def mean_intensity(self):
         return np.mean(self.intensity_image[self.image])
 
-    @_cached_property
+    @property
     def min_intensity(self):
         return np.min(self.intensity_image[self.image])
 
-    @_cached_property
+    @property
     def major_axis_length(self):
         l1, _ = self.inertia_tensor_eigvals
         return 4 * sqrt(l1)
 
-    @_cached_property
+    @property
     def minor_axis_length(self):
         _, l2 = self.inertia_tensor_eigvals
         return 4 * sqrt(l2)
@@ -243,7 +243,7 @@ def moments_central(self):
         return _moments.moments_central(self.image.astype(np.uint8),
                                         row, col, 3)
 
-    @_cached_property
+    @property
     def moments_hu(self):
         return _moments.moments_hu(self.moments_normalized)
 
@@ -251,7 +251,7 @@ def moments_hu(self):
     def moments_normalized(self):
         return _moments.moments_normalized(self.moments_central, 3)
 
-    @_cached_property
+    @property
     def orientation(self):
         a, b, b, c = self.inertia_tensor.flat
         b = -b
@@ -263,20 +263,20 @@ def orientation(self):
         else:
             return - 0.5 * atan2(2 * b, (a - c))
 
-    @_cached_property
+    @property
     def perimeter(self):
         return perimeter(self.image, 4)
 
-    @_cached_property
+    @property
     def solidity(self):
         return self.moments[0, 0] / np.sum(self.convex_image)
 
-    @_cached_property
+    @property
     def weighted_centroid(self):
         row, col = self.weighted_local_centroid
         return row + self._slice[0].start, col + self._slice[1].start
 
-    @_cached_property
+    @property
     def weighted_local_centroid(self):
         m = self.weighted_moments
         row = m[0, 1] / m[0, 0]
@@ -293,7 +293,7 @@ def weighted_moments_central(self):
         return _moments.moments_central(self._intensity_image_double,
                                         row, col, 3)
 
-    @_cached_property
+    @property
     def weighted_moments_hu(self):
         return _moments.moments_hu(self.weighted_moments_normalized)
 
diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py
index 4244738fff3..0ff5c8f727c 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 e65e756b8ad..8ab162cfb3e 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 55b279a2e54..663e1e02da9 100644
--- a/skimage/morphology/misc.py
+++ b/skimage/morphology/misc.py
@@ -15,11 +15,8 @@
 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 de1a87948c5..dbbcad8b7c2 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 d52f92bbf20..6d8097b49e4 100644
--- a/skimage/morphology/tests/test_binary.py
+++ b/skimage/morphology/tests/test_binary.py
@@ -4,7 +4,6 @@
 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 @@
 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 911e0c71005..576404f2391 100644
--- a/skimage/morphology/tests/test_grey.py
+++ b/skimage/morphology/tests/test_grey.py
@@ -5,7 +5,7 @@
 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():
-
-    def setUp(self):
-        k = 5
-        arrname = '%03i' % k
-
-        self.disk = selem.disk(k)
-
-        fname_opening = os.path.join(data_dir, "disk-open-matlab-output.npz")
-        self.expected_opening = np.load(fname_opening)[arrname]
-
-        fname_closing = os.path.join(data_dir, "disk-close-matlab-output.npz")
-        self.expected_closing = np.load(fname_closing)[arrname]
-
-    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_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)
+# 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]])
+
+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]])
+
+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]])
+
+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]])
+
+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_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_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/restoration/tests/test_unwrap.py b/skimage/restoration/tests/test_unwrap.py
index e628f4fe340..f5ff2c76930 100644
--- a/skimage/restoration/tests/test_unwrap.py
+++ b/skimage/restoration/tests/test_unwrap.py
@@ -131,6 +131,7 @@ def test_mask():
         # The end of the unwrapped array should have value equal to the
         # endpoint of the unmasked ramp
         assert_array_almost_equal_nulp(image_unwrapped[:, -1], image[i, -1])
+        assert np.ma.isMaskedArray(image_unwrapped)
 
         # Same tests, but forcing use of the 3D unwrapper by reshaping
         with expected_warnings(['length 1 dimension']):
@@ -138,7 +139,7 @@ def test_mask():
             image_wrapped_3d = image_wrapped.reshape(shape)
             image_unwrapped_3d = unwrap_phase(image_wrapped_3d)
             # remove phase shift
-            image_unwrapped_3d -= image_unwrapped_3d[0, 0, 0]  
+            image_unwrapped_3d -= image_unwrapped_3d[0, 0, 0]
         assert_array_almost_equal_nulp(image_unwrapped_3d[:, :, -1], image[i, -1])
 
 
@@ -156,5 +157,51 @@ def test_unwrap_3d_middle_wrap_around():
     unwrap = unwrap_phase(image, wrap_around=[False, True, False])
     assert np.all(unwrap == 0)
 
+
+def test_unwrap_2d_compressed_mask():
+    # ValueError when image is masked array with a compressed mask (no masked
+    # elments).  GitHub issue #1346
+    image = np.ma.zeros((10, 10))
+    unwrap = unwrap_phase(image)
+    assert np.all(unwrap == 0)
+
+
+def test_unwrap_2d_all_masked():
+    # Segmentation fault when image is masked array with a all elements masked
+    # GitHub issue #1347
+    # all elements masked
+    image = np.ma.zeros((10, 10))
+    image[:] = np.ma.masked
+    unwrap = unwrap_phase(image)
+    assert np.ma.isMaskedArray(unwrap)
+    assert np.all(unwrap.mask)
+
+    # 1 unmasked element, still zero edges
+    image = np.ma.zeros((10, 10))
+    image[:] = np.ma.masked
+    image[0, 0] = 0
+    unwrap = unwrap_phase(image)
+    assert np.ma.isMaskedArray(unwrap)
+    assert np.sum(unwrap.mask) == 99    # all but one masked
+    assert unwrap[0, 0] == 0
+
+
+def test_unwrap_3d_all_masked():
+    # all elements masked
+    image = np.ma.zeros((10, 10, 10))
+    image[:] = np.ma.masked
+    unwrap = unwrap_phase(image)
+    assert np.ma.isMaskedArray(unwrap)
+    assert np.all(unwrap.mask)
+
+    # 1 unmasked element, still zero edges
+    image = np.ma.zeros((10, 10, 10))
+    image[:] = np.ma.masked
+    image[0, 0, 0] = 0
+    unwrap = unwrap_phase(image)
+    assert np.ma.isMaskedArray(unwrap)
+    assert np.sum(unwrap.mask) == 999   # all but one masked
+    assert unwrap[0, 0, 0] == 0
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/skimage/restoration/unwrap.py b/skimage/restoration/unwrap.py
index c68102296d4..14f1eafa2f9 100644
--- a/skimage/restoration/unwrap.py
+++ b/skimage/restoration/unwrap.py
@@ -88,13 +88,14 @@ def unwrap_phase(image, wrap_around=False, seed=None):
                       'algorithm')
 
     if np.ma.isMaskedArray(image):
-        mask = np.require(image.mask, np.uint8, ['C'])
-        image = image.data
+        mask = np.require(np.ma.getmaskarray(image), np.uint8, ['C'])
     else:
         mask = np.zeros_like(image, dtype=np.uint8, order='C')
 
-    image_not_masked = np.asarray(image, dtype=np.double, order='C')
-    image_unwrapped = np.empty_like(image, dtype=np.double, order='C')
+    image_not_masked = np.asarray(
+        np.ma.getdata(image), dtype=np.double, order='C')
+    image_unwrapped = np.empty_like(image, dtype=np.double, order='C',
+                                    subok=False)
 
     if image.ndim == 1:
         unwrap_1d(image_not_masked, image_unwrapped)
diff --git a/skimage/restoration/unwrap_2d_ljmu.c b/skimage/restoration/unwrap_2d_ljmu.c
index 605c7993584..690a8d77411 100644
--- a/skimage/restoration/unwrap_2d_ljmu.c
+++ b/skimage/restoration/unwrap_2d_ljmu.c
@@ -718,10 +718,11 @@ void unwrap2D(double *wrapped_image, double *UnwrappedImage,
   horizontalEDGEs(pixel, edge, image_width, image_height, &params);
   verticalEDGEs(pixel, edge, image_width, image_height, &params);
 
-  // sort the EDGEs depending on their reiability. The PIXELs with higher
-  // relibility (small value) first
-  quicker_sort(edge, edge + params.no_of_edges - 1);
-
+  if (params.no_of_edges != 0) {
+      // sort the EDGEs depending on their reiability. The PIXELs with higher
+      // relibility (small value) first
+      quicker_sort(edge, edge + params.no_of_edges - 1);
+  }
   // gather PIXELs into groups
   gatherPIXELs(edge, &params);
 
diff --git a/skimage/restoration/unwrap_3d_ljmu.c b/skimage/restoration/unwrap_3d_ljmu.c
index be47f9260e0..fdd9847ae84 100644
--- a/skimage/restoration/unwrap_3d_ljmu.c
+++ b/skimage/restoration/unwrap_3d_ljmu.c
@@ -1136,9 +1136,11 @@ void unwrap3D(double *wrapped_volume, double *unwrapped_volume,
                 &params);
   normalEDGEs(voxel, edge, volume_width, volume_height, volume_depth, &params);
 
-  // sort the EDGEs depending on their reiability. The VOXELs with higher
-  // relibility (small value) first
-  quicker_sort(edge, edge + params.no_of_edges - 1);
+  if (params.no_of_edges != 0) {
+      // sort the EDGEs depending on their reiability. The VOXELs with higher
+      // relibility (small value) first
+      quicker_sort(edge, edge + params.no_of_edges - 1);
+  }
 
   // gather VOXELs into groups
   gatherVOXELs(edge, &params);
diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py
index 8e4ef32e1a3..1975d879b70 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 2fff52f80b5..e6e401aca6f 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 5577e46bd16..4739bd5b596 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -3,7 +3,7 @@
 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 @@
            '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 d616102d763..f8104df5f24 100644
--- a/skimage/util/arraypad.py
+++ b/skimage/util/arraypad.py
@@ -8,7 +8,7 @@
 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 eb8554b1989..d6f993f9e50 100644
--- a/skimage/util/tests/test_arraypad.py
+++ b/skimage/util/tests/test_arraypad.py
@@ -5,8 +5,8 @@
 
 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 @@ def test_check_wrong_pad_amount(self):
                       **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()
diff --git a/skimage/viewer/__init__.py b/skimage/viewer/__init__.py
index 4bf70a931e0..443e5659042 100644
--- a/skimage/viewer/__init__.py
+++ b/skimage/viewer/__init__.py
@@ -1,11 +1,6 @@
-from warnings import warn
-from skimage._shared.version_requirements import is_installed
-
+import warnings
 from .viewers import ImageViewer, CollectionViewer
-from .qt import qt_api
-
-viewer_available = not qt_api is None and is_installed('matplotlib')
-if not viewer_available:
-    warn('Viewer requires matplotlib and Qt')
+from .qt import has_qt
 
-del qt_api, is_installed, warn
+if not has_qt:
+    warnings.warn('Viewer requires Qt')
diff --git a/skimage/viewer/canvastools/base.py b/skimage/viewer/canvastools/base.py
index 197800c002e..4b1d0c00ff7 100644
--- a/skimage/viewer/canvastools/base.py
+++ b/skimage/viewer/canvastools/base.py
@@ -1,8 +1,5 @@
 import numpy as np
-try:
-    from matplotlib import lines
-except ImportError:
-    pass
+from matplotlib import lines
 
 __all__ = ['CanvasToolBase', 'ToolHandles']
 
diff --git a/skimage/viewer/canvastools/linetool.py b/skimage/viewer/canvastools/linetool.py
index f18a1915fe3..58b05bc4ca0 100644
--- a/skimage/viewer/canvastools/linetool.py
+++ b/skimage/viewer/canvastools/linetool.py
@@ -1,9 +1,6 @@
 import numpy as np
 
-try:
-    from matplotlib import lines
-except ImportError:
-    pass
+from matplotlib import lines
 from skimage.viewer.canvastools.base import CanvasToolBase, ToolHandles
 
 
diff --git a/skimage/viewer/canvastools/painttool.py b/skimage/viewer/canvastools/painttool.py
index 953f6ebe2f1..9f4bb8ecad7 100644
--- a/skimage/viewer/canvastools/painttool.py
+++ b/skimage/viewer/canvastools/painttool.py
@@ -1,11 +1,8 @@
 import numpy as np
-try:
-    import matplotlib.pyplot as plt
-    import matplotlib.colors as mcolors
-    LABELS_CMAP = mcolors.ListedColormap(['white', 'red', 'dodgerblue', 'gold',
-                                          'greenyellow', 'blueviolet'])
-except ImportError:
-    pass
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+LABELS_CMAP = mcolors.ListedColormap(['white', 'red', 'dodgerblue', 'gold',
+                                      'greenyellow', 'blueviolet'])
 from skimage.viewer.canvastools.base import CanvasToolBase
 
 
diff --git a/skimage/viewer/canvastools/recttool.py b/skimage/viewer/canvastools/recttool.py
index 99b8efeeff3..9b6c1c6b901 100644
--- a/skimage/viewer/canvastools/recttool.py
+++ b/skimage/viewer/canvastools/recttool.py
@@ -1,8 +1,4 @@
-try:
-    from matplotlib.widgets import RectangleSelector
-except ImportError:
-    RectangleSelector = object
-
+from matplotlib.widgets import RectangleSelector
 from skimage.viewer.canvastools.base import CanvasToolBase
 from skimage.viewer.canvastools.base import ToolHandles
 
diff --git a/skimage/viewer/plugins/base.py b/skimage/viewer/plugins/base.py
index def0bfd98a1..7b0418b29d1 100644
--- a/skimage/viewer/plugins/base.py
+++ b/skimage/viewer/plugins/base.py
@@ -4,19 +4,11 @@
 from warnings import warn
 
 import numpy as np
-
-from ..qt import QtGui, qt_api
-from ..qt.QtCore import Qt, Signal
+from ..qt import QtWidgets, QtCore, Signal
 from ..utils import RequiredAttr, init_qtapp
-from skimage._shared.testing import doctest_skip_parser
-
-if qt_api is not None:
-    has_qt = True
-else:
-    has_qt = False
 
 
-class Plugin(QtGui.QDialog):
+class Plugin(QtWidgets.QDialog):
     """Base class for plugins that interact with an ImageViewer.
 
     A plugin connects an image filter (or another function) to an image viewer.
@@ -101,7 +93,7 @@ def __init__(self, image_filter=None, height=0, width=400, useblit=True,
                  "then the `image_filter` argument is ignored.")
 
         self.setWindowTitle(self.name)
-        self.layout = QtGui.QGridLayout(self)
+        self.layout = QtWidgets.QGridLayout(self)
         self.resize(width, height)
         self.row = 0
 
@@ -124,7 +116,7 @@ def attach(self, image_viewer):
         the image matches the filtered value specified by attached widgets.
         """
         self.setParent(image_viewer)
-        self.setWindowFlags(Qt.Dialog)
+        self.setWindowFlags(QtCore.Qt.Dialog)
 
         self.image_viewer = image_viewer
         self.image_viewer.plugins.append(self)
diff --git a/skimage/viewer/plugins/color_histogram.py b/skimage/viewer/plugins/color_histogram.py
index 52c71786a11..d6c6ddfc2c3 100644
--- a/skimage/viewer/plugins/color_histogram.py
+++ b/skimage/viewer/plugins/color_histogram.py
@@ -1,8 +1,5 @@
 import numpy as np
-try:
-    import matplotlib.pyplot as plt
-except ImportError:
-    pass
+import matplotlib.pyplot as plt
 from skimage import color
 from skimage import exposure
 from .plotplugin import PlotPlugin
diff --git a/skimage/viewer/plugins/plotplugin.py b/skimage/viewer/plugins/plotplugin.py
index 5de6be7bb1c..f995efa8ef5 100644
--- a/skimage/viewer/plugins/plotplugin.py
+++ b/skimage/viewer/plugins/plotplugin.py
@@ -1,6 +1,6 @@
 import numpy as np
-from ..qt import QtGui
 
+from ..qt import QtGui
 from ..utils import new_plot
 from .base import Plugin
 
diff --git a/skimage/viewer/qt.py b/skimage/viewer/qt.py
new file mode 100644
index 00000000000..8f6768a51f2
--- /dev/null
+++ b/skimage/viewer/qt.py
@@ -0,0 +1,44 @@
+_qt_version = None
+has_qt = True
+
+try:
+    from matplotlib.backends.qt_compat import QtGui, QtCore, QtWidgets, QT_RC_MAJOR_VERSION as _qt_version
+except ImportError:
+    try:
+        from matplotlib.backends.qt4_compat import QtGui, QtCore
+        QtWidgets = QtGui
+        _qt_version = 4
+    except ImportError:
+        # Mock objects
+        class QtGui_cls(object):
+            QMainWindow = object
+            QDialog = object
+            QWidget = object
+
+        class QtCore_cls(object):
+            class Qt(object):
+                 TopDockWidgetArea = None
+                 BottomDockWidgetArea = None
+                 LeftDockWidgetArea = None
+                 RightDockWidgetArea = None
+
+            def Signal(self, *args, **kwargs): 
+                pass
+
+        QtGui = QtWidgets = QtGui_cls()
+        QtCore = QtCore_cls()
+
+        has_qt = False
+
+if _qt_version == 5:
+    from matplotlib.backends.backend_qt5 import FigureManagerQT
+    from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg
+elif _qt_version == 4:
+    from matplotlib.backends.backend_qt4 import FigureManagerQT
+    from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg
+else:
+    FigureManagerQT = object
+    FigureCanvasQTAgg = object
+
+Qt = QtCore.Qt
+Signal = QtCore.Signal
diff --git a/skimage/viewer/qt/QtCore.py b/skimage/viewer/qt/QtCore.py
deleted file mode 100644
index 19b92a53d1a..00000000000
--- a/skimage/viewer/qt/QtCore.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from . import qt_api
-
-if qt_api == 'pyside':
-    from PySide.QtCore import *
-elif qt_api == 'pyqt':
-    from PyQt4.QtCore import *
-    # Use pyside names for signals and slots
-    Signal = pyqtSignal
-    Slot = pyqtSlot
-else:
-    # Mock objects for buildbot (which doesn't have Qt, but imports viewer).
-    class Qt(object):
-        TopDockWidgetArea = None
-        BottomDockWidgetArea = None
-        LeftDockWidgetArea = None
-        RightDockWidgetArea = None
-
-    def Signal(*args, **kwargs):
-        pass
-
-    def Slot(*args, **kwargs):
-        pass
diff --git a/skimage/viewer/qt/QtGui.py b/skimage/viewer/qt/QtGui.py
deleted file mode 100644
index 12e9837fea9..00000000000
--- a/skimage/viewer/qt/QtGui.py
+++ /dev/null
@@ -1,11 +0,0 @@
-from . import qt_api
-
-if qt_api == 'pyside':
-    from PySide.QtGui import *
-elif qt_api == 'pyqt':
-    from PyQt4.QtGui import *
-else:
-    # Mock objects
-    QMainWindow = object
-    QDialog = object
-    QWidget = object
diff --git a/skimage/viewer/qt/README.rst b/skimage/viewer/qt/README.rst
deleted file mode 100644
index 993ae1a552c..00000000000
--- a/skimage/viewer/qt/README.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-This qt subpackage provides a wrapper to allow use of either PySide or PyQt4.
-In addition, if neither package is available, some mock objects are created to
-prevent errors in the TravisCI build. Only the objects used in the global
-namespace need to be mocked (e.g., a Qt object that gets subclassed is used
-in the global namespace).
diff --git a/skimage/viewer/qt/__init__.py b/skimage/viewer/qt/__init__.py
deleted file mode 100644
index 2fc4e9460d7..00000000000
--- a/skimage/viewer/qt/__init__.py
+++ /dev/null
@@ -1,22 +0,0 @@
-import os
-import warnings
-
-qt_api = os.environ.get('QT_API')
-
-if qt_api is None:
-    try:
-        import PyQt4
-        qt_api = 'pyqt'
-    except ImportError:
-        try:
-            import PySide
-            qt_api = 'pyside'
-        except ImportError:
-            qt_api = None
-            # Note that we don't want to raise an error because that would
-            # cause the TravisCI build to fail.
-            warnings.warn("Could not import PyQt4: ImageViewer not available!")
-
-
-if qt_api is not None:
-    os.environ['QT_API'] = qt_api
diff --git a/skimage/viewer/tests/test_plugins.py b/skimage/viewer/tests/test_plugins.py
index 528df8f4b57..8dcf7855200 100644
--- a/skimage/viewer/tests/test_plugins.py
+++ b/skimage/viewer/tests/test_plugins.py
@@ -4,14 +4,14 @@
 import skimage.data as data
 from skimage.filters.rank import median
 from skimage.morphology import disk
-from skimage.viewer import ImageViewer, viewer_available
-from numpy.testing import assert_equal, assert_allclose, assert_almost_equal
-from numpy.testing.decorators import skipif
+from skimage.viewer import ImageViewer, has_qt
+from skimage.viewer.plugins.base import Plugin
+from skimage.viewer.widgets import Slider
 from skimage.viewer.plugins import (
     LineProfile, Measure, CannyPlugin, LabelPainter, Crop, ColorHistogram,
     PlotPlugin)
-from skimage.viewer.plugins.base import Plugin
-from skimage.viewer.widgets import Slider
+from numpy.testing import assert_equal, assert_allclose, assert_almost_equal
+from numpy.testing.decorators import skipif
 from skimage._shared._warnings import expected_warnings
 
 
@@ -22,7 +22,7 @@ def setup_line_profile(image, limits='image'):
     return plugin
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_line_profile():
     """ Test a line profile using an ndim=2 image"""
     plugin = setup_line_profile(data.camera())
@@ -36,7 +36,7 @@ def test_line_profile():
     assert_allclose(scan_data.mean(), 0.2812, rtol=1e-3)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_line_profile_rgb():
     """ Test a line profile using an ndim=3 image"""
     plugin = setup_line_profile(data.chelsea(), limits=None)
@@ -51,7 +51,7 @@ def test_line_profile_rgb():
     assert_allclose(scan_data.mean(), 0.4359, rtol=1e-3)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_line_profile_dynamic():
     """Test a line profile updating after an image transform"""
     image = data.coins()[:-50, :]  # shave some off to make the line lower
@@ -77,7 +77,7 @@ def test_line_profile_dynamic():
     assert_almost_equal(np.max(line) - np.min(line), 0.639, 1)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_measure():
     image = data.camera()
     viewer = ImageViewer(image)
@@ -89,7 +89,7 @@ def test_measure():
     assert_equal(str(m._angle.text[:5]), '135.0')
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_canny():
     image = data.camera()
     viewer = ImageViewer(image)
@@ -102,7 +102,7 @@ def test_canny():
     assert edges.sum() == 2852
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_label_painter():
     image = data.camera()
     moon = data.moon()
@@ -120,7 +120,7 @@ def test_label_painter():
     assert_equal(lp.paint_tool.shape, moon.shape)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_crop():
     image = data.camera()
     viewer = ImageViewer(image)
@@ -131,7 +131,7 @@ def test_crop():
     assert_equal(viewer.image.shape, (101, 101))
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_color_histogram():
     image = skimage.img_as_float(data.load('color.png'))
     viewer = ImageViewer(image)
@@ -143,7 +143,7 @@ def test_color_histogram():
     assert_almost_equal(viewer.image.std(), 0.325, 3)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_plot_plugin():
     viewer = ImageViewer(data.moon())
     plugin = PlotPlugin(image_filter=lambda x: x)
@@ -155,7 +155,7 @@ def test_plot_plugin():
     viewer.close()
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_plugin():
     img = skimage.img_as_float(data.moon())
     viewer = ImageViewer(img)
diff --git a/skimage/viewer/tests/test_tools.py b/skimage/viewer/tests/test_tools.py
index 6b923063d2d..aae447830e7 100644
--- a/skimage/viewer/tests/test_tools.py
+++ b/skimage/viewer/tests/test_tools.py
@@ -4,14 +4,13 @@
 from numpy.testing import assert_equal
 from numpy.testing.decorators import skipif
 from skimage import data
-from skimage.viewer import ImageViewer, viewer_available
+from skimage.viewer import ImageViewer, has_qt
 from skimage.viewer.canvastools import (
     LineTool, ThickLineTool, RectangleTool, PaintTool)
 from skimage.viewer.canvastools.base import CanvasToolBase
 from matplotlib.testing.decorators import cleanup
 
 
-
 def get_end_points(image):
     h, w = image.shape[0:2]
     x = [w / 3, 2 * w / 3]
@@ -75,7 +74,7 @@ def do_event(viewer, etype, button=1, xdata=0, ydata=0, key=None):
 
 
 @cleanup
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_line_tool():
     img = data.camera()
     viewer = ImageViewer(img)
@@ -101,7 +100,7 @@ def test_line_tool():
 
 
 @cleanup
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_thick_line_tool():
     img = data.camera()
     viewer = ImageViewer(img)
@@ -125,7 +124,7 @@ def test_thick_line_tool():
 
 
 @cleanup
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_rect_tool():
     img = data.camera()
     viewer = ImageViewer(img)
@@ -154,7 +153,7 @@ def test_rect_tool():
 
 
 @cleanup
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_paint_tool():
     img = data.moon()
     viewer = ImageViewer(img)
@@ -188,7 +187,7 @@ def test_paint_tool():
 
 
 @cleanup
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_base_tool():
     img = data.moon()
     viewer = ImageViewer(img)
diff --git a/skimage/viewer/tests/test_utils.py b/skimage/viewer/tests/test_utils.py
index 77f30299017..ae158a90aa3 100644
--- a/skimage/viewer/tests/test_utils.py
+++ b/skimage/viewer/tests/test_utils.py
@@ -1,14 +1,11 @@
 # -*- coding: utf-8 -*-
-from skimage.viewer import viewer_available
-from skimage.viewer.qt import QtCore, QtGui
 from skimage.viewer import utils
 from skimage.viewer.utils import dialogs
+from skimage.viewer.qt import QtCore, QtGui, has_qt
 from numpy.testing.decorators import skipif
-from skimage.viewer import utils
-from skimage.viewer.utils import dialogs
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_event_loop():
     utils.init_qtapp()
     timer = QtCore.QTimer()
@@ -16,7 +13,7 @@ def test_event_loop():
     utils.start_qtapp()
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_format_filename():
     fname = dialogs._format_filename(('apple', 2))
     assert fname == 'apple'
@@ -24,7 +21,7 @@ def test_format_filename():
     assert fname is None
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_open_file_dialog():
     utils.init_qtapp()
     timer = QtCore.QTimer()
@@ -33,7 +30,7 @@ def test_open_file_dialog():
     assert filename is None
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_save_file_dialog():
     utils.init_qtapp()
     timer = QtCore.QTimer()
diff --git a/skimage/viewer/tests/test_viewer.py b/skimage/viewer/tests/test_viewer.py
index 1604ca6d373..548565a43ff 100644
--- a/skimage/viewer/tests/test_viewer.py
+++ b/skimage/viewer/tests/test_viewer.py
@@ -1,9 +1,11 @@
 
 from skimage import data
-from skimage.viewer.qt import QtGui, QtCore
-from skimage.viewer import ImageViewer, CollectionViewer, viewer_available
-from skimage.transform import pyramid_gaussian
+
+from skimage.viewer.qt import QtGui, QtCore, has_qt
+from skimage.viewer import ImageViewer, CollectionViewer
 from skimage.viewer.plugins import OverlayPlugin
+
+from skimage.transform import pyramid_gaussian
 from skimage.filters import sobel
 from numpy.testing import assert_equal
 from numpy.testing.decorators import skipif
@@ -11,7 +13,7 @@
 from skimage._shared._warnings import expected_warnings
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_viewer():
     astro = data.astronaut()
     coins = data.coins()
@@ -38,7 +40,7 @@ def make_key_event(key):
                            QtCore.Qt.NoModifier)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_collection_viewer():
 
     img = data.astronaut()
@@ -54,7 +56,7 @@ def test_collection_viewer():
     view._format_coord(10, 10)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 @skipif(not is_installed('matplotlib', '>=1.2'))
 def test_viewer_with_overlay():
     img = data.coins()
diff --git a/skimage/viewer/tests/test_widgets.py b/skimage/viewer/tests/test_widgets.py
index 170c186caf1..66b86fa4d76 100644
--- a/skimage/viewer/tests/test_widgets.py
+++ b/skimage/viewer/tests/test_widgets.py
@@ -1,11 +1,13 @@
 
 import os
 from skimage import data, img_as_float, io, img_as_uint
-from skimage.viewer import ImageViewer, viewer_available
+
+from skimage.viewer import ImageViewer
+from skimage.viewer.qt import QtGui, QtCore, has_qt
 from skimage.viewer.widgets import (
     Slider, OKCancelButtons, SaveButtons, ComboBox, CheckBox, Text)
 from skimage.viewer.plugins.base import Plugin
-from skimage.viewer.qt import QtGui, QtCore
+
 from numpy.testing import assert_almost_equal, assert_equal
 from numpy.testing.decorators import skipif
 from skimage._shared._warnings import expected_warnings
@@ -18,7 +20,7 @@ def get_image_viewer():
     return viewer
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_check_box():
     viewer = get_image_viewer()
     cb = CheckBox('hello', value=True, alignment='left')
@@ -33,7 +35,7 @@ def test_check_box():
     assert_equal(cb.val, False)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_combo_box():
     viewer = get_image_viewer()
     cb = ComboBox('hello', ('a', 'b', 'c'))
@@ -46,7 +48,7 @@ def test_combo_box():
     assert_equal(cb.index, 2)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_text_widget():
     viewer = get_image_viewer()
     txt = Text('hello', 'hello, world!')
@@ -57,7 +59,7 @@ def test_text_widget():
     assert_equal(str(txt.text), 'goodbye, world!')
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_slider_int():
     viewer = get_image_viewer()
     sld = Slider('radius', 2, 10, value_type='int')
@@ -71,7 +73,7 @@ def test_slider_int():
     assert_equal(sld.val, 5)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_slider_float():
     viewer = get_image_viewer()
     sld = Slider('alpha', 2.1, 3.1, value=2.1, value_type='float',
@@ -86,7 +88,7 @@ def test_slider_float():
     assert_almost_equal(sld.val, 2.5, 2)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_save_buttons():
     viewer = get_image_viewer()
     sv = SaveButtons()
@@ -114,7 +116,7 @@ def test_save_buttons():
     os.remove(filename)
 
 
-@skipif(not viewer_available)
+@skipif(not has_qt)
 def test_ok_buttons():
     viewer = get_image_viewer()
     ok = OKCancelButtons()
diff --git a/skimage/viewer/utils/core.py b/skimage/viewer/utils/core.py
index 242f2de4f48..524e9b6166c 100644
--- a/skimage/viewer/utils/core.py
+++ b/skimage/viewer/utils/core.py
@@ -1,28 +1,15 @@
 import warnings
 
 import numpy as np
-
-from ..qt import qt_api
-
-try:
-    import matplotlib as mpl
-    from matplotlib.figure import Figure
-    from matplotlib import _pylab_helpers
-    from matplotlib.colors import LinearSegmentedColormap
-    if not qt_api is None:
-        from matplotlib.backends.backend_qt4 import FigureManagerQT
-        from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg
-        if 'agg' not in mpl.get_backend().lower():
-            print("Recommended matplotlib backend is `Agg` for full "
+from skimage.viewer.qt import QtWidgets, has_qt, FigureManagerQT, FigureCanvasQTAgg
+import matplotlib as mpl
+from matplotlib.figure import Figure
+from matplotlib import _pylab_helpers
+from matplotlib.colors import LinearSegmentedColormap
+
+if has_qt and 'agg' not in mpl.get_backend().lower():
+    warnings.warn("Recommended matplotlib backend is `Agg` for full "
                   "skimage.viewer functionality.")
-    else:
-        FigureCanvasQTAgg = object
-        LinearSegmentedColormap = object
-except ImportError:
-    FigureCanvasQTAgg = object  # hack to prevent nosetest and autodoc errors
-    LinearSegmentedColormap = object
-
-from ..qt import QtGui
 
 
 __all__ = ['init_qtapp', 'start_qtapp', 'RequiredAttr', 'figimage',
@@ -39,9 +26,9 @@ def init_qtapp():
     The QApplication needs to be initialized before creating any QWidgets
     """
     global QApp
-    QApp = QtGui.QApplication.instance()
+    QApp = QtWidgets.QApplication.instance()
     if QApp is None:
-        QApp = QtGui.QApplication([])
+        QApp = QtWidgets.QApplication([])
     return QApp
 
 
@@ -128,8 +115,8 @@ def __init__(self, figure, **kwargs):
         self.fig = figure
         FigureCanvasQTAgg.__init__(self, self.fig)
         FigureCanvasQTAgg.setSizePolicy(self,
-                                        QtGui.QSizePolicy.Expanding,
-                                        QtGui.QSizePolicy.Expanding)
+                                        QtWidgets.QSizePolicy.Expanding,
+                                        QtWidgets.QSizePolicy.Expanding)
         FigureCanvasQTAgg.updateGeometry(self)
 
     def resizeEvent(self, event):
diff --git a/skimage/viewer/utils/dialogs.py b/skimage/viewer/utils/dialogs.py
index f160531db89..20dd0e16b4a 100644
--- a/skimage/viewer/utils/dialogs.py
+++ b/skimage/viewer/utils/dialogs.py
@@ -1,6 +1,6 @@
 import os
 
-from ..qt import QtGui
+from skimage.viewer.qt import QtGui
 
 
 __all__ = ['open_file_dialog', 'save_file_dialog']
diff --git a/skimage/viewer/viewers/core.py b/skimage/viewer/viewers/core.py
index 36ba6ee6272..19112dc1848 100644
--- a/skimage/viewer/viewers/core.py
+++ b/skimage/viewer/viewers/core.py
@@ -1,21 +1,15 @@
 """
 ImageViewer class for viewing and interacting with images.
 """
-from ..qt import QtGui, qt_api
-from ..qt.QtCore import Qt, Signal
-
-if qt_api is not None:
-    has_qt = True
-else:
-    has_qt = False
 
+from skimage.viewer.qt import QtWidgets, Qt, Signal
 from skimage import io, img_as_float
 from skimage.util.dtype import dtype_range
 from skimage.exposure import rescale_intensity
 import numpy as np
-from .. import utils
 from ..widgets import Slider
-from ..utils import dialogs
+from ..utils import (
+    dialogs, init_qtapp, figimage, start_qtapp, update_axes_image)
 from ..plugins.base import Plugin
 
 
@@ -152,7 +146,7 @@ def on_scroll(self, event):
             tool.on_scroll(event)
 
 
-class ImageViewer(QtGui.QMainWindow):
+class ImageViewer(QtWidgets.QMainWindow):
     """Viewer for displaying images.
 
     This viewer is a simple container object that holds a Matplotlib axes
@@ -194,7 +188,7 @@ class ImageViewer(QtGui.QMainWindow):
 
     def __init__(self, image, useblit=True):
         # Start main loop
-        utils.init_qtapp()
+        init_qtapp()
         super(ImageViewer, self).__init__()
 
         #TODO: Add ImageViewer to skimage.io window manager
@@ -202,7 +196,7 @@ def __init__(self, image, useblit=True):
         self.setAttribute(Qt.WA_DeleteOnClose)
         self.setWindowTitle("Image Viewer")
 
-        self.file_menu = QtGui.QMenu('&File', self)
+        self.file_menu = QtWidgets.QMenu('&File', self)
         self.file_menu.addAction('Open file', self.open_file,
                                  Qt.CTRL + Qt.Key_O)
         self.file_menu.addAction('Save to file', self.save_to_file,
@@ -211,7 +205,7 @@ def __init__(self, image, useblit=True):
                                  Qt.CTRL + Qt.Key_Q)
         self.menuBar().addMenu(self.file_menu)
 
-        self.main_widget = QtGui.QWidget()
+        self.main_widget = QtWidgets.QWidget()
         self.setCentralWidget(self.main_widget)
 
         if isinstance(image, Plugin):
@@ -221,7 +215,7 @@ def __init__(self, image, useblit=True):
             # When plugin is started, start
             plugin._started.connect(self._show)
 
-        self.fig, self.ax = utils.figimage(image)
+        self.fig, self.ax = figimage(image)
         self.canvas = self.fig.canvas
         self.canvas.setParent(self)
         self.ax.autoscale(enable=False)
@@ -236,7 +230,7 @@ def __init__(self, image, useblit=True):
         self._update_original_image(image)
         self.plugins = []
 
-        self.layout = QtGui.QVBoxLayout(self.main_widget)
+        self.layout = QtWidgets.QVBoxLayout(self.main_widget)
         self.layout.addWidget(self.canvas)
 
         status_bar = self.statusBar()
@@ -255,7 +249,7 @@ def __add__(self, plugin):
         if plugin.dock:
             location = self.dock_areas[plugin.dock]
             dock_location = Qt.DockWidgetArea(location)
-            dock = QtGui.QDockWidget()
+            dock = QtWidgets.QDockWidget()
             dock.setWidget(plugin)
             dock.setWindowTitle(plugin.name)
             self.addDockWidget(dock_location, dock)
@@ -341,7 +335,7 @@ def show(self, main_window=True):
         """
         self._show()
         if main_window:
-            utils.start_qtapp()
+            start_qtapp()
         return [p.output() for p in self.plugins]
 
     def redraw(self):
@@ -357,7 +351,7 @@ def image(self):
     @image.setter
     def image(self, image):
         self._img = image
-        utils.update_axes_image(self._image_plot, image)
+        update_axes_image(self._image_plot, image)
 
         # update display (otherwise image doesn't fill the canvas)
         h, w = image.shape[:2]
@@ -487,7 +481,7 @@ def update_image(self, image):
         self._update_original_image(image)
 
     def keyPressEvent(self, event):
-        if type(event) == QtGui.QKeyEvent:
+        if type(event) == QtWidgets.QKeyEvent:
             key = event.key()
             # Number keys (code: 0 = key 48, 9 = key 57) move to deciles
             if 48 <= key < 58:
diff --git a/skimage/viewer/widgets/core.py b/skimage/viewer/widgets/core.py
index 4da8925bb28..ae8be4e7ae0 100644
--- a/skimage/viewer/widgets/core.py
+++ b/skimage/viewer/widgets/core.py
@@ -15,17 +15,15 @@
         property of the same name that updates the display.
 
 """
-from ..qt import QtGui
-from ..qt import QtCore
-from ..qt.QtCore import Qt
-
+from ..qt import QtWidgets, QtCore, Qt
 from ..utils import RequiredAttr
 
 
 __all__ = ['BaseWidget', 'Slider', 'ComboBox', 'CheckBox', 'Text']
 
 
-class BaseWidget(QtGui.QWidget):
+
+class BaseWidget(QtWidgets.QWidget):
 
     plugin = RequiredAttr("Widget is not attached to a Plugin.")
 
@@ -49,11 +47,11 @@ class Text(BaseWidget):
 
     def __init__(self, name=None, text=''):
         super(Text, self).__init__(name)
-        self._label = QtGui.QLabel()
+        self._label = QtWidgets.QLabel()
         self.text = text
-        self.layout = QtGui.QHBoxLayout(self)
+        self.layout = QtWidgets.QHBoxLayout(self)
         if name is not None:
-            name_label = QtGui.QLabel()
+            name_label = QtWidgets.QLabel()
             name_label.setText(name)
             self.layout.addWidget(name_label)
         self.layout.addWidget(self._label)
@@ -105,17 +103,17 @@ def __init__(self, name, low=0.0, high=1.0, value=None, value_type='float',
         # Set widget orientation
         #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         if orientation == 'vertical':
-            self.slider = QtGui.QSlider(Qt.Vertical)
+            self.slider = QtWidgets.QSlider(Qt.Vertical)
             alignment = QtCore.Qt.AlignHCenter
             align_text = QtCore.Qt.AlignHCenter
             align_value = QtCore.Qt.AlignHCenter
-            self.layout = QtGui.QVBoxLayout(self)
+            self.layout = QtWidgets.QVBoxLayout(self)
         elif orientation == 'horizontal':
-            self.slider = QtGui.QSlider(Qt.Horizontal)
+            self.slider = QtWidgets.QSlider(Qt.Horizontal)
             alignment = QtCore.Qt.AlignVCenter
             align_text = QtCore.Qt.AlignLeft
             align_value = QtCore.Qt.AlignRight
-            self.layout = QtGui.QHBoxLayout(self)
+            self.layout = QtWidgets.QHBoxLayout(self)
         else:
             msg = "Unexpected value %s for 'orientation'"
             raise ValueError(msg % orientation)
@@ -151,11 +149,11 @@ def __init__(self, name, low=0.0, high=1.0, value=None, value_type='float',
             raise ValueError("Unexpected value %s for 'update_on'" % update_on)
         self.slider.setFocusPolicy(QtCore.Qt.StrongFocus)
 
-        self.name_label = QtGui.QLabel()
+        self.name_label = QtWidgets.QLabel()
         self.name_label.setText(self.name)
         self.name_label.setAlignment(align_text)
 
-        self.editbox = QtGui.QLineEdit()
+        self.editbox = QtWidgets.QLineEdit()
         self.editbox.setMaximumWidth(max_edit_width)
         self.editbox.setText(self.value_fmt % self.val)
         self.editbox.setAlignment(align_value)
@@ -229,20 +227,18 @@ class ComboBox(BaseWidget):
     def __init__(self, name, items, ptype='kwarg', callback=None):
         super(ComboBox, self).__init__(name, ptype, callback)
 
-        self.name_label = QtGui.QLabel()
+        self.name_label = QtWidgets.QLabel()
         self.name_label.setText(self.name)
         self.name_label.setAlignment(QtCore.Qt.AlignLeft)
 
-        self._combo_box = QtGui.QComboBox()
+        self._combo_box = QtWidgets.QComboBox()
         self._combo_box.addItems(list(items))
 
-        self.layout = QtGui.QHBoxLayout(self)
+        self.layout = QtWidgets.QHBoxLayout(self)
         self.layout.addWidget(self.name_label)
         self.layout.addWidget(self._combo_box)
 
         self._combo_box.currentIndexChanged.connect(self._value_changed)
-        # self.connect(self._combo_box,
-                # SIGNAL("currentIndexChanged(int)"), self.updateUi)
 
     @property
     def val(self):
@@ -283,11 +279,11 @@ def __init__(self, name, value=False, alignment='center', ptype='kwarg',
                  callback=None):
         super(CheckBox, self).__init__(name, ptype, callback)
 
-        self._check_box = QtGui.QCheckBox()
+        self._check_box = QtWidgets.QCheckBox()
         self._check_box.setChecked(value)
         self._check_box.setText(self.name)
 
-        self.layout = QtGui.QHBoxLayout(self)
+        self.layout = QtWidgets.QHBoxLayout(self)
         if alignment == 'center':
             self.layout.setAlignment(QtCore.Qt.AlignCenter)
         elif alignment == 'left':
diff --git a/skimage/viewer/widgets/history.py b/skimage/viewer/widgets/history.py
index d1df8689a78..ba0f0cd3402 100644
--- a/skimage/viewer/widgets/history.py
+++ b/skimage/viewer/widgets/history.py
@@ -1,8 +1,6 @@
 from textwrap import dedent
 
-from ..qt import QtGui
-from ..qt import QtCore
-
+from ..qt import QtGui, QtCore
 import numpy as np
 
 import skimage