diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
index 0907544b..48a252b2 100644
--- a/skimage/transform/__init__.py
+++ b/skimage/transform/__init__.py
@@ -4,5 +4,6 @@ from .finite_radon_transform import *
 from .integral import *
 from ._geometric import (warp, warp_coords, estimate_transform,
                          SimilarityTransform, AffineTransform,
-                         ProjectiveTransform, PolynomialTransform)
-from ._warps import resize, rotate, swirl, homography
+                         ProjectiveTransform, PolynomialTransform, 
+                         PiecewiseAffineTransform)
+from ._warps import swirl, homography
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 0f76b50e..030197da 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1,6 +1,6 @@
 import math
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage, spatial
 from skimage.util import img_as_float
 from ._warps_cy import _warp_fast
 
@@ -580,12 +580,90 @@ class PolynomialTransform(GeometricTransform):
             'parameters by exchanging source and destination coordinates,'
             'then apply the forward transformation.')
 
+class PiecewiseAffineTransform(ProjectiveTransform):
+
+    """2D piecewise affine transformation.
+
+    Control points are used to define the mapping. The transform is based on
+    a Delaunay triangulation of the points to form a mesh. Each triangle is
+    used to find a local affine transform.
+
+    Parameters
+    ----------
+    TODO
+
+    """
+
+    def __init__(self):
+        self.tess = None
+        self.triAffines = []
+        self._matrix = None
+
+    def estimate(self, src, dst):
+        """Set the control points with which to perform the piecewise affine mapping.
+
+        Number of source and destination coordinates must match.
+
+        Parameters
+        ----------
+        src : (N, 2) array
+            Source coordinates.
+        dst : (N, 2) array
+            Destination coordinates.
+
+        """
+
+        #Triangulate input positions into mesh
+        self.tess = spatial.Delaunay(src)
+
+        #Find affine mapping from source positions to destination
+        self.triAffines = []
+        for tri in self.tess.vertices:
+            affine = AffineTransform()
+            affine.estimate(src[tri,:], dst[tri,:])
+            self.triAffines.append(affine)
+
+    def __call__(self, coords):
+        """Apply forward transformation.
+
+        Parameters
+        ----------
+        coords : (N, 2) array
+            source coordinates
+
+        Returns
+        -------
+        coords : (N, 2) array
+            Transformed coordinates.
+
+        """
+        
+        out = np.ones((coords.shape[0], 2)) * -1
+
+        for ptNum, pt in enumerate(coords):
+            #Determine which triangle contains the point
+            simplexIndex = self.tess.find_simplex(pt)
+            
+            if simplexIndex == -1:
+                #This point is outside the hull of the control points
+                out[ptNum,0] = -1
+                out[ptNum,1] = -1
+                continue
+
+            #Calculate affine transformed position
+            affine = self.triAffines[simplexIndex]
+            destPos = affine(pt)
+            out[ptNum,0] = destPos[0][0]
+            out[ptNum,1] = destPos[0][1]
+
+        return out
 
 TRANSFORMS = {
     'similarity': SimilarityTransform,
     'affine': AffineTransform,
     'projective': ProjectiveTransform,
     'polynomial': PolynomialTransform,
+    'piecewiseaffine': PiecewiseAffineTransform,
 }
 HOMOGRAPHY_TRANSFORMS = (
     SimilarityTransform,
@@ -593,7 +671,6 @@ HOMOGRAPHY_TRANSFORMS = (
     ProjectiveTransform
 )
 
-
 def estimate_transform(ttype, src, dst, **kwargs):
     """Estimate 2D geometric transformation parameters.
 
@@ -604,7 +681,7 @@ def estimate_transform(ttype, src, dst, **kwargs):
 
     Parameters
     ----------
-    ttype : {'similarity', 'affine', 'projective', 'polynomial'}
+    ttype : {'similarity', 'affine', 'piecewiseaffine', 'projective', 'polynomial'}
         Type of transform.
     kwargs : array or int
         Function parameters (src, dst, n, angle)::
@@ -612,6 +689,7 @@ def estimate_transform(ttype, src, dst, **kwargs):
             NAME / TTYPE        FUNCTION PARAMETERS
             'similarity'        `src, `dst`
             'affine'            `src, `dst`
+            'piecewiseaffine'   `src, `dst`
             'projective'        `src, `dst`
             'polynomial'        `src, `dst`, `order` (polynomial order)