diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index 5d41a063..7562b650 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -684,9 +684,9 @@ def gray2rgb(image):
     """
     if np.squeeze(image).ndim == 3 and image.shape[2] in (3, 4):
         return image
-    elif is_gray(image):
+    elif image.ndim != 1 and np.squeeze(image).ndim in (1, 2, 3):
         image = image[..., np.newaxis]
-        return np.concatenate((image,)*3, axis=-1)
+        return np.concatenate(3 * (image,), axis=-1)
     else:
         raise ValueError("Input image expected to be RGB, RGBA or gray.")
 
diff --git a/skimage/filter/_denoise.py b/skimage/filter/_denoise.py
index 87850759..ff89b78b 100644
--- a/skimage/filter/_denoise.py
+++ b/skimage/filter/_denoise.py
@@ -37,7 +37,7 @@ def _denoise_tv_chambolle_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
     >>> mask = (x - 22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
     >>> mask = mask.astype(np.float)
     >>> mask += 0.2 * np.random.randn(*mask.shape)
-    >>> res = denoise_tv(mask, weight=100)
+    >>> res = denoise_tv_chambolle(mask, weight=100)
 
     """
 
@@ -127,7 +127,7 @@ def _denoise_tv_chambolle_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
     >>> from skimage import color, data
     >>> lena = color.rgb2gray(data.lena())
     >>> lena += 0.5 * lena.std() * np.random.randn(*lena.shape)
-    >>> denoised_lena = denoise_tv(lena, weight=60)
+    >>> denoised_lena = denoise_tv_chambolle(lena, weight=60)
 
     """
 
@@ -227,7 +227,7 @@ def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200,
     >>> from skimage import color, data
     >>> lena = color.rgb2gray(data.lena())
     >>> lena += 0.5 * lena.std() * np.random.randn(*lena.shape)
-    >>> denoised_lena = denoise_tv(lena, weight=60)
+    >>> denoised_lena = denoise_tv_chambolle(lena, weight=60)
 
     3D example on synthetic data:
 
@@ -235,7 +235,7 @@ def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200,
     >>> mask = (x - 22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
     >>> mask = mask.astype(np.float)
     >>> mask += 0.2*np.random.randn(*mask.shape)
-    >>> res = denoise_tv(mask, weight=100)
+    >>> res = denoise_tv_chambolle(mask, weight=100)
 
     """
 
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 23359a93..25a13765 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -951,11 +951,11 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     ----------
     image : 2-D or 3-D array
         Input image.
-    inverse_map : transformation object, callable ``xy = f(xy, **kwargs)``
+    inverse_map : transformation object, callable ``xy = f(xy, **kwargs)``, (3, 3) array
         Inverse coordinate map. A function that transforms a (N, 2) array of
         ``(x, y)`` coordinates in the *output image* into their corresponding
         coordinates in the *source image* (e.g. a transformation object or its
-        inverse).
+        inverse). See example section for usage.
     map_args : dict, optional
         Keyword arguments passed to `inverse_map`.
     output_shape : tuple (rows, cols), optional
@@ -976,26 +976,46 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
 
+    Notes
+    -----
+    In case of a `SimilarityTransform`, `AffineTransform` and
+    `ProjectiveTransform` and `order` in [0, 3] this function uses the
+    underlying transformation matrix to warp the image with a much faster
+    routine.
+
     Examples
     --------
-    Shift an image to the right:
-
     >>> from skimage.transform import warp
     >>> from skimage import data
     >>> image = data.camera()
-    >>>
-    >>> def shift_right(xy):
-    ...     xy[:, 0] -= 10
-    ...     return xy
-    >>>
-    >>> warp(image, shift_right)
 
-    Use a geometric transform to warp an image:
+    The following image warps are all equal but differ substantially in
+    execution time.
+
+    Use a geometric transform to warp an image (fast):
 
     >>> from skimage.transform import SimilarityTransform
-    >>> tform = SimilarityTransform(scale=0.1, rotation=0.1)
+    >>> tform = SimilarityTransform(translation=(0, -10))
     >>> warp(image, tform)
 
+    Shift an image to the right with a callable (slow):
+
+    >>> def shift(xy):
+    ...     xy[:, 1] -= 10
+    ...     return xy
+    >>> warp(image, shift_right)
+
+    Use a transformation matrix to warp an image (fast):
+
+    >>> matrix = np.array([[1, 0, 0], [0, 1, -10], [0, 0, 1]])
+    >>> warp(image, matrix)
+    >>> from skimage.transform import ProjectiveTransform
+    >>> warp(image, ProjectiveTransform(matrix=matrix))
+
+    You can also use the inverse of a geometric transformation (fast):
+
+    >>> warp(image, tform.inverse)
+
     """
     # Backward API compatibility
     if reverse_map is not None:
@@ -1015,16 +1035,21 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     if order in range(4) and not map_args:
         matrix = None
 
-        if inverse_map in HOMOGRAPHY_TRANSFORMS:
+        if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3):
+            matrix = inverse_map
+
+        elif inverse_map in HOMOGRAPHY_TRANSFORMS:
             matrix = inverse_map._matrix
 
-        elif hasattr(inverse_map, '__name__') \
-             and inverse_map.__name__ == 'inverse' \
-             and get_bound_method_class(inverse_map) in HOMOGRAPHY_TRANSFORMS:
+        elif (hasattr(inverse_map, '__name__')
+              and inverse_map.__name__ == 'inverse'
+              and get_bound_method_class(inverse_map)
+                  in HOMOGRAPHY_TRANSFORMS):
 
             matrix = np.linalg.inv(six.get_method_self(inverse_map)._matrix)
 
         if matrix is not None:
+            matrix = matrix.astype(np.double)
             # transform all bands
             dims = []
             for dim in range(image.shape[2]):
@@ -1042,25 +1067,30 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
         rows, cols = output_shape[:2]
 
+        if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3):
+            inverse_map = ProjectiveTransform(matrix=inverse_map)
+
         def coord_map(*args):
             return inverse_map(*args, **map_args)
 
         coords = warp_coords(coord_map, (rows, cols, bands))
 
-        # Prefilter not necessary for order 1 interpolation
+        # Prefilter not necessary for order 0, 1 interpolation
         prefilter = order > 1
         out = ndimage.map_coordinates(image, coords, prefilter=prefilter,
                                       mode=mode, order=order, cval=cval)
 
-    # The spline filters sometimes return results outside [0, 1],
-    # so clip to ensure valid data
-    clipped = np.clip(out, 0, 1)
+        # The spline filters sometimes return results outside [0, 1],
+        # so clip to ensure valid data
+        clipped = np.clip(out, 0, 1)
 
-    if mode == 'constant' and not (0 <= cval <= 1):
-        clipped[out == cval] = cval
+        if mode == 'constant' and not (0 <= cval <= 1):
+            clipped[out == cval] = cval
 
-    if clipped.ndim == 3 and orig_ndim == 2:
-        # remove singleton dim introduced by atleast_3d
-        return clipped[..., 0]
+        out = clipped
+
+    if out.ndim == 3 and orig_ndim == 2:
+        # remove singleton dimension introduced by atleast_3d
+        return out[..., 0]
     else:
-        return clipped
+        return out
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index af2b95da..7f7ef47d 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -11,10 +11,9 @@ from skimage import transform as tf, data, img_as_float
 from skimage.color import rgb2gray
 
 
-def test_warp():
-    x = np.zeros((5, 5), dtype=np.uint8)
-    x[2, 2] = 255
-    x = img_as_float(x)
+def test_warp_tform():
+    x = np.zeros((5, 5), dtype=np.double)
+    x[2, 2] = 1
     theta = - np.pi / 2
     tform = SimilarityTransform(scale=1, rotation=theta, translation=(0, 4))
 
@@ -25,10 +24,36 @@ def test_warp():
     assert_array_almost_equal(x90, np.rot90(x))
 
 
+def test_warp_callable():
+    x = np.zeros((5, 5), dtype=np.double)
+    x[2, 2] = 1
+    refx = np.zeros((5, 5), dtype=np.double)
+    refx[1, 1] = 1
+
+    shift = lambda xy: xy + 1
+
+    outx = warp(x, shift, order=1)
+    assert_array_almost_equal(outx, refx)
+
+
+def test_warp_matrix():
+    x = np.zeros((5, 5), dtype=np.double)
+    x[2, 2] = 1
+    refx = np.zeros((5, 5), dtype=np.double)
+    refx[1, 1] = 1
+
+    matrix = np.array([[1, 0, 1], [0, 1, 1], [0, 0, 1]])
+
+    # _warp_fast
+    outx = warp(x, matrix, order=1)
+    assert_array_almost_equal(outx, refx)
+    # check for ndimage.map_coordinates
+    outx = warp(x, matrix, order=5)
+
+
 def test_homography():
-    x = np.zeros((5, 5), dtype=np.uint8)
-    x[1, 1] = 255
-    x = img_as_float(x)
+    x = np.zeros((5, 5), dtype=np.double)
+    x[1, 1] = 1
     theta = -np.pi / 2
     M = np.array([[np.cos(theta), - np.sin(theta), 0],
                   [np.sin(theta),   np.cos(theta), 4],
diff --git a/skimage/util/_regular_grid.py b/skimage/util/_regular_grid.py
index e304be20..898a4aed 100644
--- a/skimage/util/_regular_grid.py
+++ b/skimage/util/_regular_grid.py
@@ -3,7 +3,7 @@ import numpy as np
 
 def regular_grid(ar_shape, n_points):
     """Find `n_points` regularly spaced along `ar_shape`.
-    
+
     The returned points (as slices) should be as close to cubically-spaced as
     possible. Essentially, the points are spaced by the Nth root of the input
     array size, where N is the number of dimensions. However, if an array
@@ -13,7 +13,7 @@ def regular_grid(ar_shape, n_points):
     Parameters
     ----------
     ar_shape : array-like of ints
-        The shape of the space embedding the grid. `len(ar_shape)` is the
+        The shape of the space embedding the grid. ``len(ar_shape)`` is the
         number of dimensions.
     n_points : int
         The (approximate) number of points to embed in the space.
@@ -66,7 +66,7 @@ def regular_grid(ar_shape, n_points):
                 break
     starts = stepsizes // 2
     stepsizes = np.round(stepsizes)
-    slices = [slice(start, None, step) for 
+    slices = [slice(start, None, step) for
               start, step in zip(starts, stepsizes)]
     slices = [slices[i] for i in unsort_dim_idxs]
     return slices
diff --git a/skimage/util/montage.py b/skimage/util/montage.py
index bdafe6ed..805b78a8 100644
--- a/skimage/util/montage.py
+++ b/skimage/util/montage.py
@@ -10,7 +10,7 @@ def montage2d(arr_in, fill='mean', rescale_intensity=False, grid_shape=None):
     """Create a 2-dimensional 'montage' from a 3-dimensional input array
     representing an ensemble of equally shaped 2-dimensional images.
 
-    For example, montage2d(arr_in, fill) with the following `arr_in`
+    For example, ``montage2d(arr_in, fill)`` with the following `arr_in`
 
     +---+---+---+
     | 1 | 2 | 3 |
@@ -37,7 +37,8 @@ def montage2d(arr_in, fill='mean', rescale_intensity=False, grid_shape=None):
     rescale_intensity: bool, optional
         Whether to rescale the intensity of each image to [0, 1].
     grid_shape: tuple, optional
-        The desired grid shape for the montage (tiles_y, tiles_x). Tthe default aspect ratio is square.
+        The desired grid shape for the montage (tiles_y, tiles_x).
+        The default aspect ratio is square.
 
     Returns
     -------
diff --git a/skimage/util/shape.py b/skimage/util/shape.py
index 1d606b42..5fe27a36 100644
--- a/skimage/util/shape.py
+++ b/skimage/util/shape.py
@@ -208,11 +208,11 @@ def view_as_windows(arr_in, window_shape, step=1):
 
     # -- basic checks on arguments
     if not isinstance(arr_in, np.ndarray):
-        raise TypeError("'arr_in' must be a numpy ndarray")
+        raise TypeError("`arr_in` must be a numpy ndarray")
     if not isinstance(window_shape, tuple):
-        raise TypeError("'window_shape' must be a tuple")
+        raise TypeError("`window_shape` must be a tuple")
     if not (len(window_shape) == arr_in.ndim):
-        raise ValueError("'window_shape' is incompatible with 'arr_in.shape'")
+        raise ValueError("`window_shape` is incompatible with `arr_in.shape`")
 
     if step < 1:
         raise ValueError("`step` must be >= 1")
@@ -221,10 +221,10 @@ def view_as_windows(arr_in, window_shape, step=1):
     window_shape = np.array(window_shape, dtype=arr_shape.dtype)
 
     if ((arr_shape - window_shape) < 0).any():
-        raise ValueError("'window_shape' is too large")
+        raise ValueError("`window_shape` is too large")
 
     if ((window_shape - 1) < 0).any():
-        raise ValueError("'window_shape' is too small")
+        raise ValueError("`window_shape` is too small")
 
     # -- build rolling window view
     arr_in = np.ascontiguousarray(arr_in)
diff --git a/skimage/util/unique.py b/skimage/util/unique.py
index e65c05ce..16a83f6b 100644
--- a/skimage/util/unique.py
+++ b/skimage/util/unique.py
@@ -9,12 +9,12 @@ def unique_rows(ar):
 
     Parameters
     ----------
-    ar : 2D np.ndarray
+    ar : 2-D ndarray
         The input array.
 
     Returns
     -------
-    ar_out : 2D np.ndarray
+    ar_out : 2-D ndarray
         A copy of the input array with repeated rows removed.
 
     Raises