diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx
index 137004e5..83722aba 100644
--- a/skimage/_shared/interpolation.pyx
+++ b/skimage/_shared/interpolation.pyx
@@ -18,8 +18,8 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols,
         Shape of image.
     r, c : int
         Position at which to interpolate.
-    mode : {'C', 'W', 'M'}
-        Wrapping mode. Constant, Wrap or Mirror.
+    mode : {'C', 'W', 'R'}
+        Wrapping mode. Constant, Wrap or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -42,8 +42,8 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
         Shape of image.
     r, c : int
         Position at which to interpolate.
-    mode : {'C', 'W', 'M'}
-        Wrapping mode. Constant, Wrap or Mirror.
+    mode : {'C', 'W', 'R'}
+        Wrapping mode. Constant, Wrap or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -76,8 +76,8 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
         Shape of image.
     r, c : int
         Position at which to get the pixel.
-    mode : {'C', 'W', 'M'}
-        Wrapping mode. Constant, Wrap or Mirror.
+    mode : {'C', 'W', 'R'}
+        Wrapping mode. Constant, Wrap or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -101,13 +101,13 @@ cdef inline int coord_map(int dim, int coord, char mode):
         Maximum coordinate.
     coord : int
         Coord provided by user.  May be < 0 or > dim.
-    mode : {'W', 'M'}
-        Whether to wrap or mirror the coordinate if it
+    mode : {'W', 'R'}
+        Whether to wrap or reflect the coordinate if it
         falls outside [0, dim).
 
     """
     dim = dim - 1
-    if mode == 'M': # mirror
+    if mode == 'R': # reflect
         if (coord < 0):
             # How many times times does the coordinate wrap?
             if (<int>(-coord / dim) % 2 != 0):
diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
index b0eee512..0735267a 100644
--- a/skimage/transform/__init__.py
+++ b/skimage/transform/__init__.py
@@ -1,9 +1,8 @@
 from .hough_transform import *
 from .radon_transform import *
 from .finite_radon_transform import *
-from ._project import homography as fast_homography
 from .integral import *
-from ._geometric import (warp, warp_coords, estimate_transform,
+from ._geometric import (estimate_transform,
                          SimilarityTransform, AffineTransform,
                          ProjectiveTransform, PolynomialTransform)
-from ._warps import swirl, homography
+from ._warps import warp, warp_coords, swirl, homography
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index c3b47a3b..e4799375 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1,30 +1,5 @@
 import math
 import numpy as np
-from scipy import ndimage
-from skimage.util import img_as_float
-
-
-def _stackcopy(a, b):
-    """Copy b into each color layer of a, such that::
-
-        a[:,:,0] = a[:,:,1] = ... = b
-
-    Parameters
-    ----------
-    a : (M, N) or (M, N, P) ndarray
-        Target array.
-    b : (M, N)
-        Source array.
-
-    Notes
-    -----
-    Color images are stored as an ``(M, N, 3)`` or ``(M, N, 4)`` arrays.
-
-    """
-    if a.ndim == 3:
-        a[:] = b[:, :, np.newaxis]
-    else:
-        a[:] = b
 
 
 class GeometricTransform(object):
@@ -603,7 +578,7 @@ class PolynomialTransform(GeometricTransform):
             'then apply the forward transformation.')
 
 
-TRANSFORMATIONS = {
+TRANSFORMS = {
     'similarity': SimilarityTransform,
     'affine': AffineTransform,
     'projective': ProjectiveTransform,
@@ -669,11 +644,11 @@ def estimate_transform(ttype, src, dst, **kwargs):
 
     """
     ttype = ttype.lower()
-    if ttype not in TRANSFORMATIONS:
+    if ttype not in TRANSFORMS:
         raise ValueError('the transformation type \'%s\' is not'
                          'implemented' % ttype)
 
-    tform = TRANSFORMATIONS[ttype]()
+    tform = TRANSFORMS[ttype]()
     tform.estimate(src, dst, **kwargs)
 
     return tform
@@ -696,158 +671,3 @@ def matrix_transform(coords, matrix):
 
     """
     return ProjectiveTransform(matrix)(coords)
-
-
-def warp_coords(orows, ocols, bands, coord_transform_fn,
-                dtype=np.float64):
-    """Build the source coordinates for the output pixels of an image warp.
-
-    Parameters
-    ----------
-    orows : int
-        Number of output rows.
-    ocols : int
-        Number of output columns.
-    bands : int
-        Number of color bands (aka channels).
-    coord_transform_fn : callable like GeometricTransform.inverse
-        Return input coordinates for given output coordinates.
-    dtype : np.dtype or string
-        dtype for return value (sane choices: float32 or float64)
-
-    Returns
-    -------
-    coords : (3, orows, ocols, bands) array of dtype `dtype`
-            Coordinates for `scipy.ndimage.map_coordinates`, that will yield
-            an image of shape (orows, ocols, bands) by drawing from source
-            points according to the `coord_transform_fn`.
-
-    Notes
-    -----
-    This is a lower-level routine that produces the source coordinates used by
-    `warp()`.
-
-    It is provided separately from `warp` to give additional flexibility to
-    users who would like, for example, to re-use a particular coordinate
-    mapping, to use specific dtypes at various points along the the
-    image-warping process, or to implement different post-processing logic
-    than `warp` performs after the call to `ndimage.map_coordinates`.
-
-
-    Examples
-    --------
-    Produce a coordinate map that Shifts an image to the right:
-
-    >>> from skimage import data
-    >>> from scipy.ndimage import map_coordinates
-    >>>
-    >>> def shift_right(xy):
-    ...     xy[:, 0] -= 10
-    ...     return xy
-    >>>
-    >>> coords = warp_coords(30, 30, 3, shift_right)
-    >>> image = data.lena().astype(np.float32)
-    >>> warped_image = map_coordinates(image, coords)
-
-    """
-
-    coords = np.empty((3, orows, ocols, bands), dtype=dtype)
-
-    # Reshape grid coordinates into a (P, 2) array of (x, y) pairs
-    tf_coords = np.indices((ocols, orows), dtype=dtype).reshape(2, -1).T
-
-    # Map each (x, y) pair to the source image according to
-    # the user-provided mapping
-    tf_coords = coord_transform_fn(tf_coords)
-
-    # Reshape back to a (2, M, N) coordinate grid
-    tf_coords = tf_coords.T.reshape((-1, ocols, orows)).swapaxes(1, 2)
-
-    # Place the y-coordinate mapping
-    _stackcopy(coords[1, ...], tf_coords[0, ...])
-
-    # Place the x-coordinate mapping
-    _stackcopy(coords[0, ...], tf_coords[1, ...])
-
-    # colour-coordinate mapping
-    coords[2, ...] = range(bands)
-
-    return coords
-
-
-def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
-         mode='constant', cval=0., reverse_map=None):
-    """Warp an image according to a given coordinate transformation.
-
-    Parameters
-    ----------
-    image : 2-D array
-        Input image.
-    inverse_map : transformation object, callable xy = f(xy, **kwargs)
-        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).
-    map_args : dict, optional
-        Keyword arguments passed to `inverse_map`.
-    output_shape : tuple (rows, cols)
-        Shape of the output image generated.
-    order : int
-        Order of splines used in interpolation. See
-        `scipy.ndimage.map_coordinates` for detail.
-    mode : string
-        How to handle values outside the image borders.  See
-        `scipy.ndimage.map_coordinates` for detail.
-    cval : float
-        Used in conjunction with mode 'constant', the value outside
-        the image boundaries.
-
-    Examples
-    --------
-    Shift an image to the right:
-
-    >>> from skimage import data
-    >>> image = data.camera()
-    >>>
-    >>> def shift_right(xy):
-    ...     xy[:, 0] -= 10
-    ...     return xy
-    >>>
-    >>> warp(image, shift_right)
-
-    """
-    # Backward API compatibility
-    if reverse_map is not None:
-        inverse_map = reverse_map
-
-    if image.ndim < 2:
-        raise ValueError("Input must have more than 1 dimension.")
-
-    image = np.atleast_3d(img_as_float(image))
-    ishape = np.array(image.shape)
-    bands = ishape[2]
-
-    if output_shape is None:
-        output_shape = ishape
-
-    rows, cols = output_shape[:2]
-
-    def coord_transform_fn(*args):
-        return inverse_map(*args, **map_args)
-
-    coords = warp_coords(rows, cols, bands, coord_transform_fn)
-
-    # Prefilter not necessary for order 1 interpolation
-    prefilter = order > 1
-    mapped = 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(mapped, 0, 1)
-
-    if mode == 'constant' and not (0 <= cval <= 1):
-        clipped[mapped == cval] = cval
-
-    # Remove singleton dim introduced by atleast_3d
-    return clipped.squeeze()
diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py
index 4c499d6b..db7e5eba 100644
--- a/skimage/transform/_warps.py
+++ b/skimage/transform/_warps.py
@@ -1,5 +1,215 @@
-from ._geometric import warp, ProjectiveTransform
 import numpy as np
+from scipy import ndimage
+from skimage.util import img_as_float
+from ._geometric import (SimilarityTransform, AffineTransform,
+                         ProjectiveTransform)
+from ._warps_cy import _warp_fast
+
+
+HOMOGRAPHY_TRANSFORMS = (
+    SimilarityTransform,
+    AffineTransform,
+    ProjectiveTransform
+)
+
+
+def _stackcopy(a, b):
+    """Copy b into each color layer of a, such that::
+
+      a[:,:,0] = a[:,:,1] = ... = b
+
+    Parameters
+    ----------
+    a : (M, N) or (M, N, P) ndarray
+        Target array.
+    b : (M, N)
+        Source array.
+
+    Notes
+    -----
+    Color images are stored as an ``(M, N, 3)`` or ``(M, N, 4)`` arrays.
+
+    """
+    if a.ndim == 3:
+        a[:] = b[:, :, np.newaxis]
+    else:
+        a[:] = b
+
+
+def warp_coords(coord_map, shape, dtype=np.float64):
+    """Build the source coordinates for the output pixels of an image warp.
+
+    Parameters
+    ----------
+    coord_map : callable like GeometricTransform.inverse
+        Return input coordinates for given output coordinates.
+    shape : tuple
+        Shape of output image ``(rows, cols[, bands])``.
+    dtype : np.dtype or string
+        dtype for return value (sane choices: float32 or float64).
+
+    Returns
+    -------
+    coords : (ndim, rows, cols[, bands]) array of dtype `dtype`
+            Coordinates for `scipy.ndimage.map_coordinates`, that will yield
+            an image of shape (orows, ocols, bands) by drawing from source
+            points according to the `coord_transform_fn`.
+
+    Notes
+    -----
+    This is a lower-level routine that produces the source coordinates used by
+    `warp()`.
+
+    It is provided separately from `warp` to give additional flexibility to
+    users who would like, for example, to re-use a particular coordinate
+    mapping, to use specific dtypes at various points along the the
+    image-warping process, or to implement different post-processing logic
+    than `warp` performs after the call to `ndimage.map_coordinates`.
+
+
+    Examples
+    --------
+    Produce a coordinate map that Shifts an image to the right:
+
+    >>> from skimage import data
+    >>> from scipy.ndimage import map_coordinates
+    >>>
+    >>> def shift_right(xy):
+    ...     xy[:, 0] -= 10
+    ...     return xy
+    >>>
+    >>> coords = warp_coords(30, 30, 3, shift_right)
+    >>> image = data.lena().astype(np.float32)
+    >>> warped_image = map_coordinates(image, coords)
+
+    """
+    rows, cols = shape[0], shape[1]
+    coords_shape = [len(shape), rows, cols]
+    if len(shape) == 3:
+        coords_shape.append(shape[2])
+    coords = np.empty(coords_shape, dtype=dtype)
+
+    # Reshape grid coordinates into a (P, 2) array of (x, y) pairs
+    tf_coords = np.indices((cols, rows), dtype=dtype).reshape(2, -1).T
+
+    # Map each (x, y) pair to the source image according to
+    # the user-provided mapping
+    tf_coords = coord_map(tf_coords)
+
+    # Reshape back to a (2, M, N) coordinate grid
+    tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
+
+    # Place the y-coordinate mapping
+    _stackcopy(coords[1, ...], tf_coords[0, ...])
+
+    # Place the x-coordinate mapping
+    _stackcopy(coords[0, ...], tf_coords[1, ...])
+
+    if len(shape) == 3:
+        coords[2, ...] = range(shape[2])
+
+    return coords
+
+
+def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
+         mode='constant', cval=0., reverse_map=None):
+    """Warp an image according to a given coordinate transformation.
+
+    Parameters
+    ----------
+    image : 2-D array
+        Input image.
+    inverse_map : transformation object, callable xy = f(xy, **kwargs)
+        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).
+    map_args : dict, optional
+        Keyword arguments passed to `inverse_map`.
+    output_shape : tuple (rows, cols)
+        Shape of the output image generated.
+    order : int
+        Order of splines used in interpolation. See
+        `scipy.ndimage.map_coordinates` for detail.
+    mode : string
+        How to handle values outside the image borders.  See
+        `scipy.ndimage.map_coordinates` for detail.
+    cval : float
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+
+    Examples
+    --------
+    Shift an image to the right:
+
+    >>> from skimage import data
+    >>> image = data.camera()
+    >>>
+    >>> def shift_right(xy):
+    ...     xy[:, 0] -= 10
+    ...     return xy
+    >>>
+    >>> warp(image, shift_right)
+
+    """
+    # Backward API compatibility
+    if reverse_map is not None:
+        inverse_map = reverse_map
+
+    if image.ndim < 2:
+        raise ValueError("Input must have more than 1 dimension.")
+
+    orig_ndim = image.ndim
+    image = np.atleast_3d(img_as_float(image))
+    ishape = np.array(image.shape)
+    bands = ishape[2]
+
+    # use fast Cython version for specific parameters
+    fast_modes = ('constant', 'reflect', 'wrap')
+    if order in (0, 1) and mode in fast_modes and not map_args:
+        matrix = None
+        if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS):
+            matrix = inverse_map._matrix
+        elif inverse_map.__name__ == 'inverse' \
+                and inverse_map.im_class in HOMOGRAPHY_TRANSFORMS:
+            matrix = np.linalg.inv(inverse_map.im_self._matrix)
+        if matrix is not None:
+            # transform all bands
+            dims = []
+            for dim in range(image.shape[2]):
+                dims.append(_warp_fast(image[..., dim], matrix,
+                            output_shape=output_shape,
+                            order=order, mode=mode, cval=cval))
+            out = np.dstack(dims)
+            if orig_ndim == 2:
+                out = out[..., 0]
+            return out
+
+    if output_shape is None:
+        output_shape = ishape
+
+    rows, cols = output_shape[:2]
+
+    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 = order > 1
+    mapped = 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(mapped, 0, 1)
+
+    if mode == 'constant' and not (0 <= cval <= 1):
+        clipped[mapped == cval] = cval
+
+    # Remove singleton dim introduced by atleast_3d
+    return clipped.squeeze()
+
 
 def _swirl_mapping(xy, center, rotation, strength, radius):
     x, y = xy.T
diff --git a/skimage/transform/_project.pyx b/skimage/transform/_warps_cy.pyx
similarity index 90%
rename from skimage/transform/_project.pyx
rename to skimage/transform/_warps_cy.pyx
index 201785ad..68ccef27 100644
--- a/skimage/transform/_project.pyx
+++ b/skimage/transform/_warps_cy.pyx
@@ -33,7 +33,7 @@ cdef inline _matrix_transform(double x, double y, double* H, double *x_,
     y_[0] = yy / zz
 
 
-def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
+def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
                mode='constant', double cval=0):
     """Projective transformation (homography).
 
@@ -71,7 +71,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
         Order of interpolation::
         * 0: Nearest-neighbour interpolation.
         * 1: Bilinear interpolation (default).
-    mode : {'constant', 'mirror', 'wrap'}
+    mode : {'constant', 'reflect', 'wrap'}
         How to handle values outside the image borders.
     cval : string
         Used in conjunction with mode 'C' (constant), the value
@@ -82,18 +82,18 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
     cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \
          np.ascontiguousarray(image, dtype=np.double)
     cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \
-         np.ascontiguousarray(np.linalg.inv(H))
+         np.ascontiguousarray(H)
 
-    if mode not in ('constant', 'wrap', 'mirror'):
+    if mode not in ('constant', 'wrap', 'reflect'):
         raise ValueError("Invalid mode specified.  Please use "
-                         "`constant`, `wrap` or `mirror`.")
+                         "`constant`, `wrap` or `reflect`.")
     cdef char mode_c
     if mode == 'constant':
         mode_c = ord('C')
     elif mode == 'wrap':
         mode_c = ord('W')
-    elif mode == 'mirror':
-        mode_c = ord('M')
+    elif mode == 'reflect':
+        mode_c = ord('R')
 
     cdef int out_r, out_c
     if output_shape is None:
@@ -116,9 +116,9 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
             _matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
             if order == 0:
                 out[tfr, tfc] = nearest_neighbour(<double*>img.data, rows,
-                                                  cols, r, c, mode_c)
+                                                  cols, r, c, mode_c, cval)
             elif order == 1:
                 out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
-                                                       cols, r, c, mode_c)
+                                                       cols, r, c, mode_c, cval)
 
     return out
diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py
index 75210b54..0e415bad 100644
--- a/skimage/transform/setup.py
+++ b/skimage/transform/setup.py
@@ -14,12 +14,12 @@ def configuration(parent_package='', top_path=None):
     config.add_data_dir('tests')
 
     cython(['_hough_transform.pyx'], working_path=base_path)
-    cython(['_project.pyx'], working_path=base_path)
+    cython(['_warps_cy.pyx'], working_path=base_path)
 
     config.add_extension('_hough_transform', sources=['_hough_transform.c'],
                          include_dirs=[get_numpy_include_dirs()])
 
-    config.add_extension('_project', sources=['_project.c'],
+    config.add_extension('_warps_cy', sources=['_warps_cy.c'],
                          include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
     return config
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index 7be3151b..f896a452 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -2,7 +2,7 @@ from numpy.testing import assert_array_almost_equal, run_module_suite
 import numpy as np
 from scipy.ndimage import map_coordinates
 
-from skimage.transform import (warp, warp_coords, fast_homography,
+from skimage.transform import (warp, warp_coords,
                                AffineTransform,
                                ProjectiveTransform,
                                SimilarityTransform)
@@ -55,11 +55,12 @@ def test_fast_homography():
     H[:2, 2] = [tx, ty]
 
     tform = ProjectiveTransform(H)
+    coords = warp_coords(tform.inverse, (img.shape[0], img.shape[1]))
 
     for order in range(2):
-        for mode in ('constant', 'mirror', 'wrap'):
-            p0 = warp(img, tform.inverse, mode=mode, order=order)
-            p1 = fast_homography(img, H, mode=mode, order=order)
+        for mode in ('constant', 'reflect', 'wrap'):
+            p0 = map_coordinates(img, coords, mode=mode, order=order)
+            p1 = warp(img, tform, mode=mode, order=order)
 
             # import matplotlib.pyplot as plt
             # f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
@@ -85,8 +86,9 @@ def test_swirl():
 
 def test_const_cval_out_of_range():
     img = np.random.randn(100, 100)
-    warped = warp(img, AffineTransform(translation=(10, 10)), cval=-10)
-    assert np.sum(warped < 0) == (2 * 100 * 10 - 10 * 10)
+    cval = - 10
+    warped = warp(img, AffineTransform(translation=(10, 10)), cval=cval)
+    assert np.sum(warped == cval) == (2 * 100 * 10 - 10 * 10)
 
 
 def test_warp_identity():
@@ -107,7 +109,7 @@ def test_warp_coords_example():
     image = data.lena().astype(np.float32)
     assert 3 == image.shape[2]
     tform = SimilarityTransform(translation=(0, -10))
-    coords = warp_coords(30, 30, 3, tform)
+    coords = warp_coords(tform, (30, 30, 3))
     warped_image1 = map_coordinates(image[:, :, 0], coords[:2])