diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index baaf064b..5961407d 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -107,3 +107,4 @@
   Polygon, circle and ellipse drawing functions
   Adaptive thresholding
   Implementation of Matlab's `regionprops`
+  Estimation of geometric transformation parameters
diff --git a/doc/examples/applications/plot_geometric.py b/doc/examples/applications/plot_geometric.py
new file mode 100644
index 00000000..337ecad7
--- /dev/null
+++ b/doc/examples/applications/plot_geometric.py
@@ -0,0 +1,132 @@
+"""
+===============================
+Using geometric transformations
+===============================
+
+In this example, we will see how to use geometric transformations in the context
+of image processing.
+"""
+
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage import transform as tf
+
+margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+
+"""
+Basics
+======
+
+Several different geometric transformation types are supported: similarity,
+affine, projective and polynomial.
+
+Geometric transformations can either be created using the explicit parameters
+(e.g. scale, shear, rotation and translation) or the transformation matrix:
+
+First we create a transformation using explicit parameters:
+"""
+
+tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 2,
+                               translation=(0, 1))
+print tform._matrix
+
+"""
+Alternatively you can define a transformation by the transformation matrix
+itself:
+"""
+
+matrix = tform._matrix.copy()
+matrix[1, 2] = 2
+tform2 = tf.SimilarityTransform(matrix)
+
+"""
+These transformation objects can then be used to apply forward and inverse
+coordinate transformations between the source and destination coordinate
+systems:
+"""
+
+coord = [1, 0]
+print tform2(coord)
+print tform2.inverse(tform(coord))
+
+"""
+Image warping
+=============
+
+Geometric transformations can also be used to warp images:
+"""
+
+text = data.text()
+
+tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 4,
+                               translation=(text.shape[0] / 2, -100))
+
+rotated = tf.warp(text, tform)
+back_rotated = tf.warp(rotated, tform.inverse)
+
+fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3))
+fig.subplots_adjust(**margins)
+plt.gray()
+ax1.imshow(text)
+ax1.axis('off')
+ax2.imshow(rotated)
+ax2.axis('off')
+ax3.imshow(back_rotated)
+ax3.axis('off')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Parameter estimation
+====================
+
+In addition to the basic functionality mentioned above you can also estimate the
+parameters of a geometric transformation using the least-squares method.
+
+This can amongst other things be used for image registration or rectification,
+where you have a set of control points or homologous/corresponding points in two
+images.
+
+Let's assume we want to recognize letters on a photograph which was not taken
+from the front but at a certain angle. In the simplest case of a plane paper
+surface the letters are projectively distorted. Simple matching algorithms would
+not be able to match such symbols. One solution to this problem would be to warp
+the image so that the distortion is removed and then apply a matching algorithm:
+"""
+
+text = data.text()
+
+src = np.array((
+    (0, 0),
+    (0, 50),
+    (300, 50),
+    (300, 0)
+))
+dst = np.array((
+    (155, 15),
+    (65, 40),
+    (260, 130),
+    (360, 95)
+))
+
+tform3 = tf.ProjectiveTransform()
+tform3.estimate(src, dst)
+warped = tf.warp(text, tform3, output_shape=(50, 300))
+
+fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(8, 3))
+fig.subplots_adjust(**margins)
+plt.gray()
+ax1.imshow(text)
+ax1.plot(dst[:, 0], dst[:, 1], '.r')
+ax1.axis('off')
+ax2.imshow(warped)
+ax2.axis('off')
+
+"""
+.. image:: PLOT2RST.current_figure
+"""
+
+plt.show()
diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
index fa4059ee..4721e7d1 100644
--- a/skimage/transform/__init__.py
+++ b/skimage/transform/__init__.py
@@ -1,8 +1,9 @@
 from .hough_transform import *
 from .radon_transform import *
 from .finite_radon_transform import *
-from .project import *
 from ._project import homography as fast_homography
 from .integral import *
-from ._warp import warp
-from ._warp_zoo import swirl
+from ._geometric import (warp, estimate_transform,
+                         SimilarityTransform, AffineTransform,
+                         ProjectiveTransform, PolynomialTransform)
+from ._warps import swirl, homography
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
new file mode 100644
index 00000000..13bcc7f6
--- /dev/null
+++ b/skimage/transform/_geometric.py
@@ -0,0 +1,774 @@
+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):
+    """Perform geometric transformations on a set of coordinates.
+
+    """
+    def __call__(self, coords):
+        """Apply forward transformation.
+
+        Parameters
+        ----------
+        coords : (N, 2) array
+            source coordinates
+
+        Returns
+        -------
+        coords : (N, 2) array
+            transformed coordinates
+
+        """
+        raise NotImplementedError()
+
+    def inverse(self, coords):
+        """Apply inverse transformation.
+
+        Parameters
+        ----------
+        coords : (N, 2) array
+            source coordinates
+
+        Returns
+        -------
+        coords : (N, 2) array
+            transformed coordinates
+
+        """
+        raise NotImplementedError()
+
+    def __add__(self, other):
+        """Combine this transformation with another.
+
+        """
+        raise NotImplementedError()
+
+
+class ProjectiveTransform(GeometricTransform):
+    """Matrix transformation.
+
+    Apply a projective transformation (homography) on coordinates.
+
+    For each homogeneous coordinate :math:`\mathbf{x} = [x, y, 1]^T`, its
+    target position is calculated by multiplying with the given matrix,
+    :math:`H`, to give :math:`H \mathbf{x}`::
+
+      [[a0 a1 a2]
+       [b0 b1 b2]
+       [c0 c1 1 ]].
+
+    E.g., to rotate by theta degrees clockwise, the matrix should be::
+
+      [[cos(theta) -sin(theta) 0]
+       [sin(theta)  cos(theta) 0]
+       [0            0         1]]
+
+    or, to translate x by 10 and y by 20::
+
+      [[1 0 10]
+       [0 1 20]
+       [0 0 1 ]].
+
+    Parameters
+    ----------
+    matrix : (3, 3) array, optional
+        Homogeneous transformation matrix.
+
+    """
+
+    coeffs = range(8)
+
+    def __init__(self, matrix=None):
+        if matrix is not None and matrix.shape != (3, 3):
+            raise ValueError("invalid shape of transformation matrix")
+        self._matrix = matrix
+
+    @property
+    def _inv_matrix(self):
+        return np.linalg.inv(self._matrix)
+
+    def _apply_mat(self, coords, matrix):
+        coords = np.array(coords, copy=False, ndmin=2)
+
+        x, y = np.transpose(coords)
+        src = np.vstack((x, y, np.ones_like(x)))
+        dst = np.dot(src.transpose(), matrix.transpose())
+
+        # rescale to homogeneous coordinates
+        dst[:, 0] /= dst[:, 2]
+        dst[:, 1] /= dst[:, 2]
+
+        return dst[:, :2]
+
+    def __call__(self, coords):
+        return self._apply_mat(coords, self._matrix)
+
+    def inverse(self, coords):
+        return self._apply_mat(coords, self._inv_matrix)
+
+    def estimate(self, src, dst):
+        """Set the transformation matrix with the explicit transformation
+        parameters.
+
+        You can determine the over-, well- and under-determined parameters
+        with the total least-squares method.
+
+        Number of source and destination coordinates must match.
+
+        The transformation is defined as::
+
+            X = (a0*x + a1*y + a2) / (c0*x + c1*y + 1)
+            Y = (b0*x + b1*y + b2) / (c0*x + c1*y + 1)
+
+        These equations can be transformed to the following form::
+
+            0 = a0*x + a1*y + a2 - c0*x*X - c1*y*X - X
+            0 = b0*x + b1*y + b2 - c0*x*Y - c1*y*Y - Y
+
+        which exist for each set of corresponding points, so we have a set of
+        N * 2 equations. The coefficients appear linearly so we can write
+        A x = 0, where::
+
+            A   = [[x y 1 0 0 0 -x*X -y*X -X]
+                   [0 0 0 x y 1 -x*Y -y*Y -Y]
+                    ...
+                    ...
+                  ]
+            x.T = [a0 a1 a2 b0 b1 b2 c0 c1 c3]
+
+        In case of total least-squares the solution of this homogeneous system
+        of equations is the right singular vector of A which corresponds to the
+        smallest singular value normed by the coefficient c3.
+
+        In case of the affine transformation the coefficients c0 and c1 are 0.
+        Thus the system of equations is::
+
+            A   = [[x y 1 0 0 0 -X]
+                   [0 0 0 x y 1 -Y]
+                    ...
+                    ...
+                  ]
+            x.T = [a0 a1 a2 b0 b1 b2 c3]
+
+        Parameters
+        ----------
+        src : (N, 2) array
+            source coordinates
+        dst : (N, 2) array
+            destination coordinates
+
+        """
+        xs = src[:, 0]
+        ys = src[:, 1]
+        xd = dst[:, 0]
+        yd = dst[:, 1]
+        rows = src.shape[0]
+
+        # params: a0, a1, a2, b0, b1, b2, c0, c1
+        A = np.zeros((rows * 2, 9))
+        A[:rows, 0] = xs
+        A[:rows, 1] = ys
+        A[:rows, 2] = 1
+        A[:rows, 6] = - xd * xs
+        A[:rows, 7] = - xd * ys
+        A[rows:, 3] = xs
+        A[rows:, 4] = ys
+        A[rows:, 5] = 1
+        A[rows:, 6] = - yd * xs
+        A[rows:, 7] = - yd * ys
+        A[:rows, 8] = xd
+        A[rows:, 8] = yd
+
+        # Select relevant columns, depending on params
+        A = A[:, self.coeffs + [8]]
+
+        _, _, V = np.linalg.svd(A)
+
+        H = np.zeros((3, 3))
+        # solution is right singular vector that corresponds to smallest
+        # singular value
+        H.flat[self.coeffs + [8]] = - V[-1, :-1] / V[-1, -1]
+        H[2, 2] = 1
+
+        self._matrix = H
+
+    def __add__(self, other):
+        """Combine this transformation with another.
+
+        """
+        if isinstance(other, ProjectiveTransform):
+            # combination of the same types result in a transformation of this
+            # type again, otherwise use general projective transformation
+            if type(self) == type(other):
+                tform = self.__class__
+            else:
+                tform = ProjectiveTransform
+            return tform(other._matrix.dot(self._matrix))
+        else:
+            raise TypeError("Cannot combine transformations of differing "
+                            "types.")
+
+
+class AffineTransform(ProjectiveTransform):
+
+    """2D affine transformation of the form::
+
+        X = a0*x + a1*y + a2 =
+          = sx*x*cos(rotation) - sy*y*sin(rotation + shear) + a2
+
+        Y = b0*x + b1*y + b2 =
+          = sx*x*sin(rotation) + sy*y*cos(rotation + shear) + b2
+
+    where ``sx`` and ``sy`` are zoom factors in the x and y directions,
+    and the homogeneous transformation matrix is::
+
+        [[a0  a1  a2]
+         [b0  b1  b2]
+         [0   0    1]]
+
+    Parameters
+    ----------
+    matrix : (3, 3) array, optional
+        Homogeneous transformation matrix.
+    scale : (sx, sy) as array, list or tuple, optional
+        scale factors
+    rotation : float, optional
+        rotation angle in counter-clockwise direction as radians
+    shear : float, optional
+        shear angle in counter-clockwise direction as radians
+    translation : (tx, ty) as array, list or tuple, optional
+        translation parameters
+
+    """
+
+    coeffs = range(6)
+
+    def __init__(self, matrix=None, scale=None, rotation=None, shear=None,
+                 translation=None):
+        self._matrix = matrix
+        params = any(param is not None
+                     for param in (scale, rotation, shear, translation))
+
+        if params and matrix is not None:
+            raise ValueError("You cannot specify the transformation matrix and "
+                             "the implicit parameters at the same time.")
+        elif matrix is not None and matrix.shape != (3, 3):
+            raise ValueError("Invalid shape of transformation matrix.")
+        elif params:
+            if scale is None:
+                scale = (1, 1)
+            if rotation is None:
+                rotation = 0
+            if shear is None:
+                shear = 0
+            if translation is None:
+                translation = (0, 0)
+
+            sx, sy = scale
+            self._matrix = np.array([
+                [sx * math.cos(rotation), - sy * math.sin(rotation + shear), 0],
+                [sx * math.sin(rotation),   sy * math.cos(rotation + shear), 0],
+                [                      0,                                 0, 1]
+            ])
+            self._matrix[0:2, 2] = translation
+
+    @property
+    def scale(self):
+        sx = math.sqrt(self._matrix[0, 0] ** 2 + self._matrix[1, 0] ** 2)
+        sy = math.sqrt(self._matrix[0, 1] ** 2 + self._matrix[1, 1] ** 2)
+        return sx, sy
+
+    @property
+    def rotation(self):
+        return math.atan2(self._matrix[1, 0], self._matrix[0, 0])
+
+    @property
+    def shear(self):
+        beta = math.atan2(- self._matrix[0, 1], self._matrix[1, 1])
+        return beta - self.rotation
+
+    @property
+    def translation(self):
+        return self._matrix[0:2, 2]
+
+
+class SimilarityTransform(ProjectiveTransform):
+    """2D similarity transformation of the form::
+
+        X = a0*x - b0*y + a1 =
+          = m*x*cos(rotation) + m*y*sin(rotation) + a1
+
+        Y = b0*x + a0*y + b1 =
+          = m*x*sin(rotation) + m*y*cos(rotation) + b1
+
+    where ``m`` is a zoom factor and the homogeneous transformation matrix is::
+
+        [[a0  b0  a1]
+         [b0  a0  b1]
+         [0   0    1]]
+
+    Parameters
+    ----------
+    matrix : (3, 3) array, optional
+        Homogeneous transformation matrix.
+    scale : float, optional
+        scale factor
+    rotation : float, optional
+        rotation angle in counter-clockwise direction as radians
+    translation : (tx, ty) as array, list or tuple, optional
+        x, y translation parameters
+
+    """
+
+    def __init__(self, matrix=None, scale=None, rotation=None,
+                 translation=None):
+        self._matrix = matrix
+        params = any(param is not None
+                     for param in (scale, rotation, translation))
+
+        if params and matrix is not None:
+            raise ValueError("You cannot specify the transformation matrix and "
+                             "the implicit parameters at the same time.")
+        elif matrix is not None and matrix.shape != (3, 3):
+            raise ValueError("Invalid shape of transformation matrix.")
+        elif params:
+            if scale is None:
+                scale = 1
+            if rotation is None:
+                rotation = 0
+            if translation is None:
+                translation = (0, 0)
+
+            self._matrix = np.array([
+                [math.cos(rotation), - math.sin(rotation), 0],
+                [math.sin(rotation),   math.cos(rotation), 0],
+                [                 0,                    0, 1]
+            ])
+            self._matrix *= scale
+            self._matrix[0:2, 2] = translation
+
+    def estimate(self, src, dst):
+        """Set the transformation matrix with the explicit parameters.
+
+        You can determine the over-, well- and under-determined parameters
+        with the total least-squares method.
+
+        Number of source and destination coordinates must match.
+
+        The transformation is defined as::
+
+            X = a0*x - b0*y + a1
+            Y = b0*x + a0*y + b1
+
+        These equations can be transformed to the following form::
+
+            0 = a0*x - b0*y + a1 - X
+            0 = b0*x + a0*y + b1 - Y
+
+        which exist for each set of corresponding points, so we have a set of
+        N * 2 equations. The coefficients appear linearly so we can write
+        A x = 0, where::
+
+            A   = [[x 1 -y 0 -X]
+                   [y 0  x 1 -Y]
+                    ...
+                    ...
+                  ]
+            x.T = [a0 a1 b0 b1 c3]
+
+        In case of total least-squares the solution of this homogeneous system
+        of equations is the right singular vector of A which corresponds to the
+        smallest singular value normed by the coefficient c3.
+
+        Parameters
+        ----------
+        src : (N, 2) array
+            source coordinates
+        dst : (N, 2) array
+            destination coordinates
+
+        """
+        xs = src[:, 0]
+        ys = src[:, 1]
+        xd = dst[:, 0]
+        yd = dst[:, 1]
+        rows = src.shape[0]
+
+        # params: a0, a1, b0, b1
+        A = np.zeros((rows * 2, 5))
+        A[:rows, 0] = xs
+        A[:rows, 2] = - ys
+        A[:rows, 1] = 1
+        A[rows:, 2] = xs
+        A[rows:, 0] = ys
+        A[rows:, 3] = 1
+        A[:rows, 4] = xd
+        A[rows:, 4] = yd
+
+        _, _, V = np.linalg.svd(A)
+
+        # solution is right singular vector that corresponds to smallest
+        # singular value
+        a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1]
+
+        self._matrix = np.array([[a0, -b0, a1],
+                                 [b0,  a0, b1],
+                                 [ 0,   0,  1]])
+
+    @property
+    def scale(self):
+        if math.cos(self.rotation) == 0:
+            # sin(self.rotation) == 1
+            scale = self._matrix[0, 1]
+        else:
+            scale = self._matrix[0, 0] / math.cos(self.rotation)
+        return scale
+
+    @property
+    def rotation(self):
+        return math.atan2(self._matrix[1, 0], self._matrix[1, 1])
+
+    @property
+    def translation(self):
+        return self._matrix[0:2, 2]
+
+
+class PolynomialTransform(GeometricTransform):
+    """2D transformation of the form::
+
+        X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i ))
+        Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i ))
+
+    Parameters
+    ----------
+    params : (2, N) array, optional
+        Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So,
+        a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`.
+
+    """
+
+    def __init__(self, params=None):
+        if params is not None and params.shape[0] != 2:
+            raise ValueError("invalid shape of transformation parameters")
+        self._params = params
+
+    def estimate(self, src, dst, order):
+        """Set the transformation matrix with the explicit transformation
+        parameters.
+
+        You can determine the over-, well- and under-determined parameters
+        with the total least-squares method.
+
+        Number of source and destination coordinates must match.
+
+        The transformation is defined as::
+
+            X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i ))
+            Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i ))
+
+        These equations can be transformed to the following form::
+
+            0 = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) - X
+            0 = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) - Y
+
+        which exist for each set of corresponding points, so we have a set of
+        N * 2 equations. The coefficients appear linearly so we can write
+        A x = 0, where::
+
+            A   = [[1 x y x**2 x*y y**2 ... 0 ...             0 -X]
+                   [0 ...                 0 1 x y x**2 x*y y**2 -Y]
+                    ...
+                    ...
+                  ]
+            x.T = [a00 a10 a11 a20 a21 a22 ... ann
+                   b00 b10 b11 b20 b21 b22 ... bnn c3]
+
+        In case of total least-squares the solution of this homogeneous system
+        of equations is the right singular vector of A which corresponds to the
+        smallest singular value normed by the coefficient c3.
+
+        Parameters
+        ----------
+        src : (N, 2) array
+            source coordinates
+        dst : (N, 2) array
+            destination coordinates
+        order : int
+            polynomial order (number of coefficients is order + 1)
+
+        """
+        xs = src[:, 0]
+        ys = src[:, 1]
+        xd = dst[:, 0]
+        yd = dst[:, 1]
+        rows = src.shape[0]
+
+        # number of unknown polynomial coefficients
+        u = (order + 1) * (order + 2)
+
+        A = np.zeros((rows * 2, u + 1))
+        pidx = 0
+        for j in range(order + 1):
+            for i in range(j + 1):
+                A[:rows, pidx] = xs ** (j - i) * ys ** i
+                A[rows:, pidx + u / 2] = xs ** (j - i) * ys ** i
+                pidx += 1
+
+        A[:rows, -1] = xd
+        A[rows:, -1] = yd
+
+        _, _, V = np.linalg.svd(A)
+
+        # solution is right singular vector that corresponds to smallest
+        # singular value
+        params = - V[-1, :-1] / V[-1, -1]
+
+        self._params = params.reshape((2, u / 2))
+
+    def __call__(self, coords):
+        """Apply forward transformation.
+
+        Parameters
+        ----------
+        coords : (N, 2) array
+            source coordinates
+
+        Returns
+        -------
+        coords : (N, 2) array
+            transformed coordinates
+
+        """
+        x = coords[:, 0]
+        y = coords[:, 1]
+        u = len(self._params.ravel())
+        # number of coefficients -> u = (order + 1) * (order + 2)
+        order = int((- 3 + math.sqrt(9 - 4 * (2 - u))) / 2)
+        dst = np.zeros(coords.shape)
+
+        pidx = 0
+        for j in range(order + 1):
+            for i in range(j + 1):
+                dst[:, 0] += self._params[0, pidx] * x ** (j - i) * y ** i
+                dst[:, 1] += self._params[1, pidx] * x ** (j - i) * y ** i
+                pidx += 1
+
+        return dst
+
+    def inverse(self, coords):
+        raise Exception(
+            'There is no explicit way to do the inverse polynomial '
+            'transformation. Instead, estimate the inverse transformation '
+            'parameters by exchanging source and destination coordinates,'
+            'then apply the forward transformation.')
+
+
+TRANSFORMATIONS = {
+    'similarity': SimilarityTransform,
+    'affine': AffineTransform,
+    'projective': ProjectiveTransform,
+    'polynomial': PolynomialTransform,
+}
+
+
+def estimate_transform(ttype, src, dst, **kwargs):
+    """Estimate 2D geometric transformation parameters.
+
+    You can determine the over-, well- and under-determined parameters
+    with the total least-squares method.
+
+    Number of source and destination coordinates must match.
+
+    Parameters
+    ----------
+    ttype : {'similarity', 'affine', 'projective', 'polynomial'}
+        Type of transform.
+    kwargs : array or int
+        Function parameters (src, dst, n, angle)::
+
+            NAME / TTYPE        FUNCTION PARAMETERS
+            'similarity'        `src, `dst`
+            'affine'            `src, `dst`
+            'projective'        `src, `dst`
+            'polynomial'        `src, `dst`, `order` (polynomial order)
+
+        Also see examples below.
+
+    Returns
+    -------
+    tform : :class:`GeometricTransform`
+        Transform object containing the transformation parameters and providing
+        access to forward and inverse transformation functions.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from skimage import transform as tf
+
+    >>> # estimate transformation parameters
+    >>> src = np.array([0, 0, 10, 10]).reshape((2, 2))
+    >>> dst = np.array([12, 14, 1, -20]).reshape((2, 2))
+
+    >>> tform = tf.estimate_transform('similarity', src, dst)
+
+    >>> tform.inverse(tform(src)) # == src
+
+    >>> # warp image using the estimated transformation
+    >>> from skimage import data
+    >>> image = data.camera()
+
+    >>> 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))
+
+    >>> # unite transformations, applied in order from left to right
+    >>> tform3 = tform + tform2
+    >>> tform3(src) # == tform2(tform(src))
+
+    """
+    ttype = ttype.lower()
+    if ttype not in TRANSFORMATIONS:
+        raise ValueError('the transformation type \'%s\' is not'
+                         'implemented' % ttype)
+
+    tform = TRANSFORMATIONS[ttype]()
+    tform.estimate(src, dst, **kwargs)
+
+    return tform
+
+
+def matrix_transform(coords, matrix):
+    """Apply 2D matrix transform.
+
+    Parameters
+    ----------
+    coords : (N, 2) array
+        x, y coordinates to transform
+    matrix : (3, 3) array
+        Homogeneous transformation matrix.
+
+    Returns
+    -------
+    coords : (N, 2) array
+            transformed coordinates
+
+    """
+    return ProjectiveTransform(matrix)(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 : string
+        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
+
+    coords = np.empty(np.r_[3, output_shape], dtype=float)
+
+    ## Construct transformed coordinates
+
+    rows, cols = output_shape[:2]
+
+    # Reshape grid coordinates into a (P, 2) array of (x, y) pairs
+    tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T
+
+    # Map each (x, y) pair to the source image according to
+    # the user-provided mapping
+    tf_coords = inverse_map(tf_coords, **map_args)
+
+    # 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, ...])
+
+    # colour-coordinate mapping
+    coords[2, ...] = range(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
+    return np.clip(mapped.squeeze(), 0, 1)
diff --git a/skimage/transform/_warp.py b/skimage/transform/_warp.py
deleted file mode 100644
index 6477c988..00000000
--- a/skimage/transform/_warp.py
+++ /dev/null
@@ -1,113 +0,0 @@
-__all__ = ['warp']
-
-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 ``MxNx3`` or ``MxNx4`` arrays.
-
-    """
-    if a.ndim == 3:
-        a[:] = b[:, :, np.newaxis]
-    else:
-        a[:] = b
-
-
-def warp(image, reverse_map, map_args={},
-         output_shape=None, order=1, mode='constant', cval=0.):
-    """Warp an image according to a given coordinate transformation.
-
-    Parameters
-    ----------
-    image : 2-D array
-        Input image.
-    reverse_map : callable xy = f(xy, **kwargs)
-        Reverse coordinate map.  A function that transforms a Px2 array of
-        ``(x, y)`` coordinates in the *output image* into their corresponding
-        coordinates in the *source image*.  Also see examples below.
-    map_args : dict, optional
-        Keyword arguments passed to `reverse_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 : string
-        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)
-
-    """
-    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
-
-    coords = np.empty(np.r_[3, output_shape], dtype=float)
-
-    ## Construct transformed coordinates
-
-    rows, cols = output_shape[:2]
-
-    # Reshape grid coordinates into a (P, 2) array of (x, y) pairs
-    tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T
-
-    # Map each (x, y) pair to the source image according to
-    # the user-provided mapping
-    tf_coords = reverse_map(tf_coords, **map_args)
-
-    # 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, ...])
-
-    # colour-coordinate mapping
-    coords[2, ...] = range(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
-    return np.clip(mapped.squeeze(), 0, 1)
diff --git a/skimage/transform/_warp_zoo.py b/skimage/transform/_warp_zoo.py
deleted file mode 100644
index 4df531d1..00000000
--- a/skimage/transform/_warp_zoo.py
+++ /dev/null
@@ -1,75 +0,0 @@
-from __future__ import division
-import numpy as np
-
-from ._warp import warp
-
-
-def _swirl_mapping(xy, center, rotation, strength, radius):
-    x, y = xy.T
-    x0, y0 = center
-    rho = np.sqrt((x - x0)**2 + (y - y0)**2)
-
-    # Ensure that the transformation decays to approximately 1/1000-th
-    # within the specified radius.
-    radius = radius / 5 * np.log(2)
-
-    theta = rotation + strength * \
-            np.exp(-rho / radius) + \
-            np.arctan2(y - y0, x - x0)
-
-    xy[..., 0] = x0 + rho * np.cos(theta)
-    xy[..., 1] = y0 + rho * np.sin(theta)
-
-    return xy
-
-
-def swirl(image, center=None, strength=1, radius=100, rotation=0,
-          output_shape=None, order=1, mode='constant', cval=0):
-    """Perform a swirl transformation.
-
-    Parameters
-    ----------
-    image : ndarray
-        Input image.
-    center : (x,y) tuple or (2,) ndarray
-        Center coordinate of transformation.
-    strength : float
-        The amount of swirling applied.
-    radius : float
-        The extent of the swirl in pixels.  The effect dies out
-        rapidly beyond `radius`.
-    rotation : float
-        Additional rotation applied to the image.
-
-    Returns
-    -------
-    swirled : ndarray
-        Swirled version of the input.
-
-    Other parameters
-    ----------------
-    output_shape : tuple or ndarray
-        Size of the generated output image.
-    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.
-
-    """
-
-    if center is None:
-        center = np.array(image.shape)[:2] / 2
-
-    warp_args = {'center': center,
-                 'rotation': rotation,
-                 'strength': strength,
-                 'radius': radius}
-
-    return warp(image, _swirl_mapping, map_args=warp_args,
-                output_shape=output_shape,
-                order=order, mode=mode, cval=cval)
diff --git a/skimage/transform/project.py b/skimage/transform/_warps.py
similarity index 50%
rename from skimage/transform/project.py
rename to skimage/transform/_warps.py
index a140e03e..f09f7944 100644
--- a/skimage/transform/project.py
+++ b/skimage/transform/_warps.py
@@ -1,14 +1,75 @@
-"""Image projection.
-
-"""
-
+from ._geometric import warp, ProjectiveTransform
 import numpy as np
-from scipy.ndimage import interpolation as ndii
-from ._warp import _stackcopy
 
-__all__ = ['homography']
+def _swirl_mapping(xy, center, rotation, strength, radius):
+    x, y = xy.T
+    x0, y0 = center
+    rho = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
 
-eps = np.finfo(float).eps
+    # Ensure that the transformation decays to approximately 1/1000-th
+    # within the specified radius.
+    radius = radius / 5 * np.log(2)
+
+    theta = rotation + strength * \
+            np.exp(-rho / radius) + \
+            np.arctan2(y - y0, x - x0)
+
+    xy[..., 0] = x0 + rho * np.cos(theta)
+    xy[..., 1] = y0 + rho * np.sin(theta)
+
+    return xy
+
+
+def swirl(image, center=None, strength=1, radius=100, rotation=0,
+          output_shape=None, order=1, mode='constant', cval=0):
+    """Perform a swirl transformation.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    center : (x,y) tuple or (2,) ndarray
+        Center coordinate of transformation.
+    strength : float
+        The amount of swirling applied.
+    radius : float
+        The extent of the swirl in pixels.  The effect dies out
+        rapidly beyond `radius`.
+    rotation : float
+        Additional rotation applied to the image.
+
+    Returns
+    -------
+    swirled : ndarray
+        Swirled version of the input.
+
+    Other parameters
+    ----------------
+    output_shape : tuple or ndarray
+        Size of the generated output image.
+    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.
+
+    """
+
+    if center is None:
+        center = np.array(image.shape)[:2] / 2
+
+    warp_args = {'center': center,
+                 'rotation': rotation,
+                 'strength': strength,
+                 'radius': radius}
+
+    return warp(image, _swirl_mapping, map_args=warp_args,
+                output_shape=output_shape,
+                order=order, mode=mode, cval=cval)
 
 
 def homography(image, H, output_shape=None, order=1,
@@ -89,49 +150,11 @@ def homography(image, H, output_shape=None, order=1,
            [  0,   0,   0,   0,   0]], dtype=uint8)
 
     """
-    if image.ndim < 2:
-        raise ValueError("Input must have more than 1 dimension.")
+    import warnings
+    warnings.warn('the homography function is deprecated; '
+                  'use the `warp` and `ProjectiveTransform` class instead',
+                  category=DeprecationWarning)
 
-    image = np.atleast_3d(image)
-    ishape = np.array(image.shape)
-    bands = ishape[2]
-
-    if output_shape is None:
-        output_shape = ishape
-
-    coords = np.empty(np.r_[3, output_shape], dtype=float)
-
-    # TODO: Refactor this method to use transform.warp instead.
-
-    # Construct transformed coordinates
-    rows, cols = output_shape[:2]
-    rows, cols = np.mgrid[:rows, :cols]
-    tf_coords = np.empty(shape=cols.shape,
-                         dtype=[('cols', float),
-                                ('rows', float),
-                                ('z', float)])
-    tf_coords['cols'], tf_coords['rows'] = cols, rows
-    tf_coords['z'] = 1
-    tf_coords = tf_coords.view((float, 3))
-
-    tf_coords = np.dot(tf_coords, np.linalg.inv(H).transpose())
-    tf_coords[np.absolute(tf_coords) < eps] = 0.
-
-    # normalize coordinates
-    tf_coords[..., :2] /= tf_coords[..., 2, np.newaxis]
-
-    # y-coordinate mapping
-    _stackcopy(coords[0, ...], tf_coords[..., 1])
-
-    # x-coordinate mapping
-    _stackcopy(coords[1, ...], tf_coords[..., 0])
-
-    # colour-coordinate mapping
-    coords[2, ...] = range(bands)
-
-    # Prefilter not necessary for order 1 interpolation
-    prefilter = order > 1
-    mapped = ndii.map_coordinates(image, coords, prefilter=prefilter,
-                                  mode=mode, order=order, cval=cval)
-
-    return mapped.squeeze()
+    tform = ProjectiveTransform(H)
+    return warp(image, inverse_map=tform.inverse, output_shape=output_shape,
+                order=order, mode=mode, cval=cval)
diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py
new file mode 100644
index 00000000..d3e7f20c
--- /dev/null
+++ b/skimage/transform/tests/test_geometric.py
@@ -0,0 +1,159 @@
+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)
+
+
+SRC = np.array([
+    [-12.3705, -10.5075],
+    [-10.7865, 15.4305],
+    [8.6985, 10.8675],
+    [11.4975, -9.5715],
+    [7.8435, 7.4835],
+    [-5.3325, 6.5025],
+    [6.7905, -6.3765],
+    [-6.1695, -0.8235],
+])
+DST = np.array([
+    [0, 0],
+    [0, 5800],
+    [4900, 5800],
+    [4900, 0],
+    [4479, 4580],
+    [1176, 3660],
+    [3754, 790],
+    [1024, 1931],
+])
+
+
+def test_stackcopy():
+    layers = 4
+    x = np.empty((3, 3, layers))
+    y = np.eye(3, 3)
+    _stackcopy(x, y)
+    for i in range(layers):
+        assert_array_almost_equal(x[..., i], y)
+
+
+def test_similarity_estimation():
+    # exact solution
+    tform = estimate_transform('similarity', SRC[:2, :], DST[:2, :])
+    assert_array_almost_equal(tform(SRC[:2, :]), DST[:2, :])
+    assert_equal(tform._matrix[0, 0], tform._matrix[1, 1])
+    assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0])
+
+    # over-determined
+    tform2 = estimate_transform('similarity', SRC, DST)
+    assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC)
+    assert_equal(tform2._matrix[0, 0], tform2._matrix[1, 1])
+    assert_equal(tform2._matrix[0, 1], - tform2._matrix[1, 0])
+
+    # via estimate method
+    tform3 = SimilarityTransform()
+    tform3.estimate(SRC, DST)
+    assert_array_almost_equal(tform3._matrix, tform2._matrix)
+
+
+def test_similarity_init():
+    # init with implicit parameters
+    scale = 0.1
+    rotation = 1
+    translation = (1, 1)
+    tform = SimilarityTransform(scale=scale, rotation=rotation,
+                                translation=translation)
+    assert_array_almost_equal(tform.scale, scale)
+    assert_array_almost_equal(tform.rotation, rotation)
+    assert_array_almost_equal(tform.translation, translation)
+
+    # init with transformation matrix
+    tform2 = SimilarityTransform(tform._matrix)
+    assert_array_almost_equal(tform2.scale, scale)
+    assert_array_almost_equal(tform2.rotation, rotation)
+    assert_array_almost_equal(tform2.translation, translation)
+
+
+def test_affine_estimation():
+    # exact solution
+    tform = estimate_transform('affine', SRC[:3, :], DST[:3, :])
+    assert_array_almost_equal(tform(SRC[:3, :]), DST[:3, :])
+
+    # over-determined
+    tform2 = estimate_transform('affine', SRC, DST)
+    assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC)
+
+    # via estimate method
+    tform3 = AffineTransform()
+    tform3.estimate(SRC, DST)
+    assert_array_almost_equal(tform3._matrix, tform2._matrix)
+
+
+def test_affine_init():
+    # init with implicit parameters
+    scale = (0.1, 0.13)
+    rotation = 1
+    shear = 0.1
+    translation = (1, 1)
+    tform = AffineTransform(scale=scale, rotation=rotation, shear=shear,
+                            translation=translation)
+    assert_array_almost_equal(tform.scale, scale)
+    assert_array_almost_equal(tform.rotation, rotation)
+    assert_array_almost_equal(tform.shear, shear)
+    assert_array_almost_equal(tform.translation, translation)
+
+    # init with transformation matrix
+    tform2 = AffineTransform(tform._matrix)
+    assert_array_almost_equal(tform2.scale, scale)
+    assert_array_almost_equal(tform2.rotation, rotation)
+    assert_array_almost_equal(tform2.shear, shear)
+    assert_array_almost_equal(tform2.translation, translation)
+
+
+def test_projective_estimation():
+    # exact solution
+    tform = estimate_transform('projective', SRC[:4, :], DST[:4, :])
+    assert_array_almost_equal(tform(SRC[:4, :]), DST[:4, :])
+
+    # over-determined
+    tform2 = estimate_transform('projective', SRC, DST)
+    assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC)
+
+    # via estimate method
+    tform3 = ProjectiveTransform()
+    tform3.estimate(SRC, DST)
+    assert_array_almost_equal(tform3._matrix, tform2._matrix)
+
+
+def test_projective_init():
+    tform = estimate_transform('projective', SRC, DST)
+    # init with transformation matrix
+    tform2 = ProjectiveTransform(tform._matrix)
+    assert_array_almost_equal(tform2._matrix, tform._matrix)
+
+
+def test_polynomial_estimation():
+    # over-determined
+    tform = estimate_transform('polynomial', SRC, DST, order=10)
+    assert_array_almost_equal(tform(SRC), DST, 6)
+
+    # via estimate method
+    tform2 = PolynomialTransform()
+    tform2.estimate(SRC, DST, order=10)
+    assert_array_almost_equal(tform2._params, tform._params)
+
+
+def test_union():
+    tform1 = SimilarityTransform(scale=0.1, rotation=0.3)
+    tform2 = SimilarityTransform(scale=0.1, rotation=0.9)
+    tform3 = SimilarityTransform(scale=0.1 ** 2, rotation=0.3 + 0.9)
+
+    tform = tform1 + tform2
+
+    assert_array_almost_equal(tform._matrix, tform3._matrix)
+
+
+if __name__ == "__main__":
+    from numpy.testing import run_module_suite
+    run_module_suite()
diff --git a/skimage/transform/tests/test_project.py b/skimage/transform/tests/test_project.py
deleted file mode 100644
index 25ebda20..00000000
--- a/skimage/transform/tests/test_project.py
+++ /dev/null
@@ -1,63 +0,0 @@
-import numpy as np
-from numpy.testing import assert_array_almost_equal
-
-from skimage.transform._warp import _stackcopy
-from skimage.transform import homography, fast_homography
-from skimage import data
-from skimage.color import rgb2gray
-
-
-def test_stackcopy():
-    layers = 4
-    x = np.empty((3, 3, layers))
-    y = np.eye(3, 3)
-    _stackcopy(x, y)
-    for i in range(layers):
-        assert_array_almost_equal(x[..., i], y)
-
-
-def test_homography():
-    x = np.arange(9, dtype=np.uint8).reshape((3, 3)) + 1
-    theta = -np.pi / 2
-    M = np.array([[np.cos(theta), -np.sin(theta), 0],
-                  [np.sin(theta),  np.cos(theta), 2],
-                  [0,              0,             1]])
-    x90 = homography(x, M, order=1)
-    assert_array_almost_equal(x90, np.rot90(x))
-
-
-def test_fast_homography():
-    img = rgb2gray(data.lena()).astype(np.uint8)
-    img = img[:, :100]
-
-    theta = np.deg2rad(30)
-    scale = 0.5
-    tx, ty = 50, 50
-
-    H = np.eye(3)
-    S = scale * np.sin(theta)
-    C = scale * np.cos(theta)
-
-    H[:2, :2] = [[C, -S], [S, C]]
-    H[:2, 2] = [tx, ty]
-
-    for mode in ('constant', 'mirror', 'wrap'):
-        p0 = homography(img, H, mode=mode, order=1)
-        p1 = fast_homography(img, H, mode=mode)
-        p1 = np.round(p1)
-
-        ## import matplotlib.pyplot as plt
-        ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
-        ## ax0.imshow(img)
-        ## ax1.imshow(p0, cmap=plt.cm.gray)
-        ## ax2.imshow(p1, cmap=plt.cm.gray)
-        ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
-        ## plt.show()
-
-        d = np.mean(np.abs(p0 - p1))
-        assert d < 0.2
-
-
-if __name__ == "__main__":
-    from numpy.testing import run_module_suite
-    run_module_suite()
diff --git a/skimage/transform/tests/test_swirl.py b/skimage/transform/tests/test_swirl.py
deleted file mode 100644
index d71f8231..00000000
--- a/skimage/transform/tests/test_swirl.py
+++ /dev/null
@@ -1,17 +0,0 @@
-import numpy as np
-from numpy.testing import assert_array_almost_equal
-
-from skimage import transform as tf, data, img_as_float
-
-
-def test_roundtrip():
-    image = img_as_float(data.checkerboard())
-
-    swirl_params = {'radius': 80, 'rotation': 0, 'order': 2, 'mode': 'reflect'}
-    swirled = tf.swirl(image, strength=10, **swirl_params)
-    unswirled = tf.swirl(swirled, strength=-10, **swirl_params)
-
-    assert np.mean(np.abs(image - unswirled)) < 0.01
-
-if __name__ == "__main__":
-    np.testing.run_module_suite()
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
new file mode 100644
index 00000000..7c2e52f2
--- /dev/null
+++ b/skimage/transform/tests/test_warps.py
@@ -0,0 +1,79 @@
+from numpy.testing import assert_array_almost_equal, run_module_suite
+import numpy as np
+
+from skimage.transform import (warp, homography, fast_homography,
+                               SimilarityTransform)
+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)
+    theta = - np.pi / 2
+    tform = SimilarityTransform(scale=1, rotation=theta, translation=(0, 4))
+
+    x90 = warp(x, tform, order=1)
+    assert_array_almost_equal(x90, np.rot90(x))
+
+    x90 = warp(x, tform.inverse, order=1)
+    assert_array_almost_equal(x90, np.rot90(x))
+
+
+def test_homography():
+    x = np.zeros((5, 5), dtype=np.uint8)
+    x[1, 1] = 255
+    x = img_as_float(x)
+    theta = -np.pi/2
+    M = np.array([[np.cos(theta),-np.sin(theta),0],
+                  [np.sin(theta), np.cos(theta),4],
+                  [0,             0,            1]])
+    x90 = homography(x, M, order=1)
+    assert_array_almost_equal(x90, np.rot90(x))
+
+
+def test_fast_homography():
+    img = rgb2gray(data.lena()).astype(np.uint8)
+    img = img[:, :100]
+
+    theta = np.deg2rad(30)
+    scale = 0.5
+    tx, ty = 50, 50
+
+    H = np.eye(3)
+    S = scale * np.sin(theta)
+    C = scale * np.cos(theta)
+
+    H[:2, :2] = [[C, -S], [S, C]]
+    H[:2, 2] = [tx, ty]
+
+    for mode in ('constant', 'mirror', 'wrap'):
+        p0 = homography(img, H, mode=mode, order=1)
+        p1 = fast_homography(img, H, mode=mode)
+        p1 = np.round(p1)
+
+        ## import matplotlib.pyplot as plt
+        ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
+        ## ax0.imshow(img)
+        ## ax1.imshow(p0, cmap=plt.cm.gray)
+        ## ax2.imshow(p1, cmap=plt.cm.gray)
+        ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
+        ## plt.show()
+
+        d = np.mean(np.abs(p0 - p1))
+        assert d < 0.2
+
+
+def test_swirl():
+    image = img_as_float(data.checkerboard())
+
+    swirl_params = {'radius': 80, 'rotation': 0, 'order': 2, 'mode': 'reflect'}
+    swirled = tf.swirl(image, strength=10, **swirl_params)
+    unswirled = tf.swirl(swirled, strength=-10, **swirl_params)
+
+    assert np.mean(np.abs(image - unswirled)) < 0.01
+
+
+if __name__ == "__main__":
+    run_module_suite()