diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd
index c883d00d..6038a853 100644
--- a/skimage/_shared/interpolation.pxd
+++ b/skimage/_shared/interpolation.pxd
@@ -1,13 +1,24 @@
 
-cdef inline double nearest_neighbour(double* image, int rows, int cols,
-                                     double r, double c, char mode,
-                                     double cval=*)
+cdef inline double nearest_neighbour_interpolation(double* image, int rows,
+                                                   int cols, double r,
+                                                   double c, char mode,
+                                                   double cval) nogil
 
 cdef inline double bilinear_interpolation(double* image, int rows, int cols,
                                           double r, double c, char mode,
-                                          double cval=*)
+                                          double cval) nogil
+
+cdef inline double quadratic_interpolation(double x, double[3] f) nogil
+cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
+                                             double r, double c, char mode,
+                                             double cval) nogil
+
+cdef inline double cubic_interpolation(double x, double[4] f) nogil
+cdef inline double bicubic_interpolation(double* image, int rows, int cols,
+                                         double r, double c, char mode,
+                                         double cval) nogil
 
 cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
-                             char mode, double cval=*)
+                             char mode, double cval) nogil
 
-cdef inline int coord_map(int dim, int coord, char mode)
+cdef inline int coord_map(int dim, int coord, char mode) nogil
diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx
index 137004e5..231ca045 100644
--- a/skimage/_shared/interpolation.pyx
+++ b/skimage/_shared/interpolation.pyx
@@ -5,24 +5,30 @@
 from libc.math cimport ceil, floor, round
 
 
-cdef inline double nearest_neighbour(double* image, int rows, int cols,
-                                     double r, double c, char mode,
-                                     double cval=0):
+cdef inline double nearest_neighbour_interpolation(double* image, int rows,
+                                                   int cols, double r,
+                                                   double c, char mode,
+                                                   double cval) nogil:
     """Nearest neighbour interpolation at a given position in the image.
 
     Parameters
     ----------
     image : double array
         Input image.
-    rows, cols: int
+    rows, cols : int
         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', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
     cval : double
         Constant value to use for constant mode.
 
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
     """
 
     return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
@@ -31,22 +37,27 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols,
 
 cdef inline double bilinear_interpolation(double* image, int rows, int cols,
                                           double r, double c, char mode,
-                                          double cval=0):
+                                          double cval) nogil:
     """Bilinear interpolation at a given position in the image.
 
     Parameters
     ----------
     image : double array
         Input image.
-    rows, cols: int
+    rows, cols : int
         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', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
     cval : double
         Constant value to use for constant mode.
 
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
     """
     cdef double dr, dc
     cdef int minr, minc, maxr, maxc
@@ -64,23 +75,172 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
     return (1 - dr) * top + dr * bottom
 
 
+cdef inline double quadratic_interpolation(double x, double[3] f) nogil:
+    """Quadratic interpolation.
+
+    Parameters
+    ----------
+    x : double
+        Position in the interval [-1, 1].
+    f : double[4]
+        Function values at positions [-1, 0, 1].
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+    return f[1] - 0.25 * (f[0] - f[2]) * x
+
+
+cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
+                                             double r, double c, char mode,
+                                             double cval) nogil:
+    """Biquadratic interpolation at a given position in the image.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols : int
+        Shape of image.
+    r, c : int
+        Position at which to interpolate.
+    mode : {'C', 'W', 'R', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    cval : double
+        Constant value to use for constant mode.
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+
+    cdef int r0 = <int>round(r)
+    cdef int c0 = <int>round(c)
+    if r < 0:
+        r0 -= 1
+    if c < 0:
+        c0 -= 1
+    # scale position to range [-1, 1]
+    cdef double xr = (r - r0) - 1
+    cdef double xc = (c - c0) - 1
+    if r == r0:
+        xr += 1
+    if c == c0:
+        xc += 1
+
+    cdef double fc[3], fr[3]
+
+    cdef int pr, pc
+
+    # row-wise cubic interpolation
+    for pr in range(r0, r0 + 3):
+        for pc in range(c0, c0 + 3):
+            fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
+        fr[pr - r0] = quadratic_interpolation(xc, fc)
+
+    # cubic interpolation for interpolated values of each row
+    return quadratic_interpolation(xr, fr)
+
+
+cdef inline double cubic_interpolation(double x, double[4] f) nogil:
+    """Cubic interpolation.
+
+    Parameters
+    ----------
+    x : double
+        Position in the interval [0, 1].
+    f : double[4]
+        Function values at positions [0, 1/3, 2/3, 1].
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+    return \
+        f[1] + 0.5 * x * \
+            (f[2] - f[0] + x * \
+                (2.0 * f[0] - 5.0 * f[1] + 4.0 * f[2] - f[3] + x * \
+                    (3.0 * (f[1] - f[2]) + f[3] - f[0])))
+
+
+cdef inline double bicubic_interpolation(double* image, int rows, int cols,
+                                         double r, double c, char mode,
+                                         double cval) nogil:
+    """Bicubic interpolation at a given position in the image.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols : int
+        Shape of image.
+    r, c : int
+        Position at which to interpolate.
+    mode : {'C', 'W', 'R', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    cval : double
+        Constant value to use for constant mode.
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+
+    cdef int r0 = <int>r - 1
+    cdef int c0 = <int>c - 1
+    if r < 0:
+        r0 -= 1
+    if c < 0:
+        c0 -= 1
+    # scale position to range [0, 1]
+    cdef double xr = (r - r0) / 3
+    cdef double xc = (c - c0) / 3
+
+    cdef double fc[4], fr[4]
+
+    cdef int pr, pc
+
+    # row-wise cubic interpolation
+    for pr in range(r0, r0 + 4):
+        for pc in range(c0, c0 + 4):
+            fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
+        fr[pr - r0] = cubic_interpolation(xc, fc)
+
+    # cubic interpolation for interpolated values of each row
+    return cubic_interpolation(xr, fr)
+
+
 cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
-                             char mode, double cval=0):
+                             char mode, double cval) nogil:
     """Get a pixel from the image, taking wrapping mode into consideration.
 
     Parameters
     ----------
     image : double array
         Input image.
-    rows, cols: int
+    rows, cols : int
         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', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
     cval : double
         Constant value to use for constant mode.
 
+    Returns
+    -------
+    value : double
+        Pixel value at given position.
+
     """
     if mode == 'C':
         if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
@@ -91,7 +251,7 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
         return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
 
 
-cdef inline int coord_map(int dim, int coord, char mode):
+cdef inline int coord_map(int dim, int coord, char mode) nogil:
     """
     Wrap a coordinate, according to a given mode.
 
@@ -101,28 +261,33 @@ 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', 'N'}
+        Whether to wrap or reflect the coordinate if it
         falls outside [0, dim).
 
     """
     dim = dim - 1
-    if mode == 'M': # mirror
-        if (coord < 0):
+    if mode == 'R': # reflect
+        if coord < 0:
             # How many times times does the coordinate wrap?
-            if (<int>(-coord / dim) % 2 != 0):
+            if <int>(-coord / dim) % 2 != 0:
                 return dim - <int>(-coord % dim)
             else:
                 return <int>(-coord % dim)
-        elif (coord > dim):
-            if (<int>(coord / dim) % 2 != 0):
+        elif coord > dim:
+            if <int>(coord / dim) % 2 != 0:
                 return <int>(dim - (coord % dim))
             else:
                 return <int>(coord % dim)
     elif mode == 'W': # wrap
-        if (coord < 0):
+        if coord < 0:
             return <int>(dim - (-coord % dim))
-        elif (coord > dim):
+        elif coord > dim:
             return <int>(coord % dim)
+    elif mode == 'N': # nearest
+        if coord < 0:
+            return 0
+        elif coord > dim:
+            return dim
 
     return coord
diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx
index 5fb53e90..70a446bb 100644
--- a/skimage/feature/_texture.pyx
+++ b/skimage/feature/_texture.pyx
@@ -132,7 +132,7 @@ def _local_binary_pattern(np.ndarray[double, ndim=2] image,
         for c in range(image.shape[1]):
             for i in range(P):
                 texture[i] = bilinear_interpolation(<double*>image.data,
-                    rows, cols, r + coords[i, 0], c + coords[i, 1], 'C')
+                    rows, cols, r + coords[i, 0], c + coords[i, 1], 'C', 0)
             # signed / thresholded texture
             for i in range(P):
                 if texture[i] - image[r, c] >= 0:
diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
index b0eee512..0907544b 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,
                          SimilarityTransform, AffineTransform,
                          ProjectiveTransform, PolynomialTransform)
-from ._warps import swirl, homography
+from ._warps import resize, rotate, swirl, homography
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index c3b47a3b..270198aa 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -2,29 +2,7 @@ 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
+from ._warps_cy import _warp_fast
 
 
 class GeometricTransform(object):
@@ -603,12 +581,17 @@ class PolynomialTransform(GeometricTransform):
             'then apply the forward transformation.')
 
 
-TRANSFORMATIONS = {
+TRANSFORMS = {
     'similarity': SimilarityTransform,
     'affine': AffineTransform,
     'projective': ProjectiveTransform,
     'polynomial': PolynomialTransform,
 }
+HOMOGRAPHY_TRANSFORMS = (
+    SimilarityTransform,
+    AffineTransform,
+    ProjectiveTransform
+)
 
 
 def estimate_transform(ttype, src, dst, **kwargs):
@@ -660,8 +643,8 @@ def estimate_transform(ttype, src, dst, **kwargs):
     >>> warp(image, inverse_map=tform.inverse)
 
     >>> # create transformation with explicit parameters
-    >>> tform2 = tf.SimilarityTransform()
-    >>> tform2.compose_implicit(scale=1.1, rotation=1, translation=(10, 20))
+    >>> tform2 = tf.SimilarityTransform(scale=1.1, rotation=1,
+    ...     translation=(10, 20))
 
     >>> # unite transformations, applied in order from left to right
     >>> tform3 = tform + tform2
@@ -669,11 +652,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
@@ -698,26 +681,44 @@ def matrix_transform(coords, matrix):
     return ProjectiveTransform(matrix)(coords)
 
 
-def warp_coords(orows, ocols, bands, coord_transform_fn,
-                dtype=np.float64):
+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
     ----------
-    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
+    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)
+        dtype for return value (sane choices: float32 or float64).
 
     Returns
     -------
-    coords : (3, orows, ocols, bands) array of dtype `dtype`
+    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`.
@@ -745,23 +746,26 @@ def warp_coords(orows, ocols, bands, coord_transform_fn,
     ...     xy[:, 0] -= 10
     ...     return xy
     >>>
-    >>> coords = warp_coords(30, 30, 3, shift_right)
     >>> image = data.lena().astype(np.float32)
+    >>> coords = warp_coords(shift_right, image.shape)
     >>> warped_image = map_coordinates(image, coords)
 
     """
-
-    coords = np.empty((3, orows, ocols, bands), dtype=dtype)
+    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((ocols, orows), dtype=dtype).reshape(2, -1).T
+    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_transform_fn(tf_coords)
+    tf_coords = coord_map(tf_coords)
 
     # Reshape back to a (2, M, N) coordinate grid
-    tf_coords = tf_coords.T.reshape((-1, ocols, orows)).swapaxes(1, 2)
+    tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
 
     # Place the y-coordinate mapping
     _stackcopy(coords[1, ...], tf_coords[0, ...])
@@ -769,8 +773,8 @@ def warp_coords(orows, ocols, bands, coord_transform_fn,
     # Place the x-coordinate mapping
     _stackcopy(coords[0, ...], tf_coords[1, ...])
 
-    # colour-coordinate mapping
-    coords[2, ...] = range(bands)
+    if len(shape) == 3:
+        coords[2, ...] = range(shape[2])
 
     return coords
 
@@ -823,31 +827,53 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     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]
 
-    if output_shape is None:
-        output_shape = ishape
+    # use fast Cython version for specific interpolation orders
+    if order in range(4) 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]
 
-    rows, cols = output_shape[:2]
+    else: # use ndimage.map_coordinates
 
-    def coord_transform_fn(*args):
-        return inverse_map(*args, **map_args)
+        if output_shape is None:
+            output_shape = ishape
 
-    coords = warp_coords(rows, cols, bands, coord_transform_fn)
+        rows, cols = output_shape[:2]
 
-    # Prefilter not necessary for order 1 interpolation
-    prefilter = order > 1
-    mapped = ndimage.map_coordinates(image, coords, prefilter=prefilter,
-                                     mode=mode, order=order, cval=cval)
+        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
+        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(mapped, 0, 1)
+    clipped = np.clip(out, 0, 1)
 
     if mode == 'constant' and not (0 <= cval <= 1):
-        clipped[mapped == cval] = cval
+        clipped[out == 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 f09f7944..1ed0bed9 100644
--- a/skimage/transform/_warps.py
+++ b/skimage/transform/_warps.py
@@ -1,5 +1,119 @@
-from ._geometric import warp, ProjectiveTransform
 import numpy as np
+from ._geometric import warp, SimilarityTransform, AffineTransform
+
+
+def resize(image, output_shape, order=1, mode='constant', cval=0.):
+    """Resize image.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    output_shape : tuple or ndarray
+        Size of the generated output image `(rows, cols)`.
+
+    Returns
+    -------
+    resized : ndarray
+        Resized version of the input.
+
+    Other parameters
+    ----------------
+    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 : string
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+
+    """
+
+    rows, cols = output_shape
+    orig_rows, orig_cols = image.shape[0], image.shape[1]
+
+    rscale = float(orig_rows) / rows
+    cscale = float(orig_cols) / cols
+
+    # 3 control points necessary to estimate exact AffineTransform
+    src_corners = np.array([[1, 1], [1, rows], [cols, rows]]) - 1
+    dst_corners = np.zeros(src_corners.shape, dtype=np.double)
+    # take into account that 0th pixel is at position (0.5, 0.5)
+    dst_corners[:, 0] = cscale * (src_corners[:, 0] + 0.5) - 0.5
+    dst_corners[:, 1] = rscale * (src_corners[:, 1] + 0.5) - 0.5
+
+    tform = AffineTransform()
+    tform.estimate(src_corners, dst_corners)
+
+    return warp(image, tform, output_shape=output_shape, order=order,
+                mode=mode, cval=cval)
+
+
+def rotate(image, angle, resize=False, order=1, mode='constant', cval=0.):
+    """Rotate image by a certain angle around its center.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    angle : float
+        Rotation angle in degrees in counter-clockwise direction.
+    resize: bool, optional
+        Determine whether the shape of the output image will be automatically
+        calculated, so the complete rotated image exactly fits. Default is
+        False.
+
+    Returns
+    -------
+    rotated : ndarray
+        Rotated version of the input.
+
+    Other parameters
+    ----------------
+    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 : string
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+
+    """
+
+    rows, cols = image.shape[0], image.shape[1]
+
+    # rotation around center
+    translation = np.array((cols, rows)) / 2. - 0.5
+    tform1 = SimilarityTransform(translation=-translation)
+    tform2 = SimilarityTransform(rotation=np.deg2rad(angle))
+    tform3 = SimilarityTransform(translation=translation)
+    tform = tform1 + tform2 + tform3
+
+    output_shape = None
+    if not resize:
+        # determine shape of output image
+        corners = np.array([[1, 1], [1, rows], [cols, rows], [cols, 1]])
+        corners = tform(corners - 1)
+        minc = corners[:, 0].min()
+        minr = corners[:, 1].min()
+        maxc = corners[:, 0].max()
+        maxr = corners[:, 1].max()
+        out_rows = maxr - minr + 1
+        out_cols = maxc - minc + 1
+        output_shape = np.ceil((out_rows, out_cols))
+
+        # fit output image in new shape
+        translation = ((cols - out_cols) / 2., (rows - out_rows) / 2.)
+        tform4 = SimilarityTransform(translation=translation)
+        tform = tform4 + tform
+
+    return warp(image, tform, output_shape=output_shape, order=order,
+                mode=mode, cval=cval)
+
 
 def _swirl_mapping(xy, center, rotation, strength, radius):
     x, y = xy.T
@@ -74,7 +188,14 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
 
 def homography(image, H, output_shape=None, order=1,
                mode='constant', cval=0.):
-    """Perform a projective transformation (homography) on an image.
+    """
+    .. note:: Deprecated in skimage 0.7
+        `homography` will be removed in skimage 0.8, it is replaced by
+        `warp` because the latter provides the same functionality::
+
+            warp(image, ProjectiveTransform(H))
+
+    Perform a projective transformation (homography) on an image.
 
     For each pixel, given its homogeneous coordinate :math:`\mathbf{x}
     = [x, y, 1]^T`, its target position is calculated by multiplying
diff --git a/skimage/transform/_project.pyx b/skimage/transform/_warps_cy.pyx
similarity index 61%
rename from skimage/transform/_project.pyx
rename to skimage/transform/_warps_cy.pyx
index 4df6e001..9d3dca70 100644
--- a/skimage/transform/_project.pyx
+++ b/skimage/transform/_warps_cy.pyx
@@ -5,12 +5,15 @@
 
 cimport numpy as np
 import numpy as np
-from skimage._shared.interpolation cimport (nearest_neighbour,
-                                            bilinear_interpolation)
+from cython.parallel import prange
+from skimage._shared.interpolation cimport (nearest_neighbour_interpolation,
+                                            bilinear_interpolation,
+                                            biquadratic_interpolation,
+                                            bicubic_interpolation)
 
 
-cdef inline _matrix_transform(double x, double y, double* H, double *x_,
-                              double *y_):
+cdef inline void _matrix_transform(double x, double y, double* H, double *x_,
+                                   double *y_) nogil:
     """Apply a homography to a coordinate.
 
     Parameters
@@ -33,10 +36,9 @@ 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).
+    """Projective transformation (homography).
 
     Perform a projective transformation (homography) of a
     floating point image, using bi-linear interpolation.
@@ -72,7 +74,9 @@ 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'}
+        * 2: Biquadratic interpolation (default).
+        * 3: Bicubic interpolation.
+    mode : {'constant', 'reflect', 'wrap'}
         How to handle values outside the image borders.
     cval : string
         Used in conjunction with mode 'C' (constant), the value
@@ -83,18 +87,12 @@ 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', 'nearest'):
         raise ValueError("Invalid mode specified.  Please use "
-                         "`constant`, `wrap` or `mirror`.")
-    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')
+                         "`constant`, `nearest`, `wrap` or `reflect`.")
+    cdef char mode_c = ord(mode[0].upper())
 
     cdef int out_r, out_c
     if output_shape is None:
@@ -104,22 +102,32 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
         out_r = output_shape[0]
         out_c = output_shape[1]
 
-    cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
+    cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] out = \
          np.zeros((out_r, out_c), dtype=np.double)
+    cdef double* out_data = <double*>out.data
 
     cdef int tfr, tfc
     cdef double r, c
     cdef int rows = img.shape[0]
     cdef int cols = img.shape[1]
 
-    for tfr in range(out_r):
+    cdef double (*interp_func)(double*, int, int, double, double,
+                               char, double) nogil
+    if order == 0:
+        interp_func = nearest_neighbour_interpolation
+    elif order == 1:
+        interp_func = bilinear_interpolation
+    elif order == 2:
+        interp_func = biquadratic_interpolation
+    elif order == 3:
+        interp_func = bicubic_interpolation
+
+    for tfr in prange(out_r, nogil=True):
+        # make r, c thread local variables
+        r = c = 0
         for tfc in range(out_c):
             _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)
-            elif order == 1:
-                out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
-                                                       cols, r, c, mode_c)
+            out_data[tfr * out_r + tfc] = interp_func(<double*>img.data, rows,
+                                                      cols, r, c, mode_c, cval)
 
     return out
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index f80cc3cb..0213b2bc 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -15,7 +15,7 @@ References:
 from __future__ import division
 import numpy as np
 from scipy.fftpack import fftshift, fft, ifft
-from ._project import homography
+from ._warps_cy import _warp_fast
 
 __all__ = ["radon", "iradon"]
 
@@ -77,8 +77,8 @@ def radon(image, theta=None):
         return shift1.dot(R).dot(shift0)
 
     for i in range(len(theta)):
-        rotated = homography(padded_image,
-                             build_rotation(-theta[i]))
+        rotated = _warp_fast(padded_image,
+                             np.linalg.inv(build_rotation(-theta[i])))
 
         out[:, i] = rotated.sum(0)[::-1]
 
diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py
index 75210b54..cd2a0f8a 100644
--- a/skimage/transform/setup.py
+++ b/skimage/transform/setup.py
@@ -14,13 +14,15 @@ 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'],
-                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
+    config.add_extension('_warps_cy', sources=['_warps_cy.c'],
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'],
+                         extra_compile_args=['-fopenmp'],
+                         extra_link_args=['-fopenmp'])
 
     return config
 
diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py
index abe8b2e7..57e95d46 100644
--- a/skimage/transform/tests/test_geometric.py
+++ b/skimage/transform/tests/test_geometric.py
@@ -1,10 +1,9 @@
 import numpy as np
 from numpy.testing import assert_equal, assert_array_almost_equal
-
 from skimage.transform._geometric import _stackcopy
-from skimage.transform import (estimate_transform, SimilarityTransform,
-                               AffineTransform, ProjectiveTransform,
-                               PolynomialTransform)
+from skimage.transform import (estimate_transform,
+                               SimilarityTransform, AffineTransform,
+                               ProjectiveTransform, PolynomialTransform)
 
 
 SRC = np.array([
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index 7be3151b..ac9272c6 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, rotate, resize,
                                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 order in range(4):
+        for mode in ('constant', 'reflect', 'wrap', 'nearest'):
+            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)
@@ -73,6 +74,22 @@ def test_fast_homography():
             assert d < 0.001
 
 
+def test_rotate():
+    x = np.zeros((5, 5), dtype=np.double)
+    x[1, 1] = 1
+    x90 = rotate(x, 90)
+    assert_array_almost_equal(x90, np.rot90(x))
+
+
+def test_resize():
+    x = np.zeros((5, 5), dtype=np.double)
+    x[1, 1] = 1
+    resized = resize(x, (10, 10), order=0)
+    ref = np.zeros((10, 10))
+    ref[2:4, 2:4] = 1
+    assert_array_almost_equal(resized, ref)
+
+
 def test_swirl():
     image = img_as_float(data.checkerboard())
 
@@ -85,8 +102,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 +125,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])