diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index d3fd10ad..f1b611eb 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -951,7 +951,7 @@ 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
@@ -1022,16 +1022,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]):
@@ -1049,12 +1054,15 @@ 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)
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],