diff --git a/doc/examples/plot_swirl.py b/doc/examples/plot_swirl.py
index c49dbabf..18947dcb 100644
--- a/doc/examples/plot_swirl.py
+++ b/doc/examples/plot_swirl.py
@@ -36,7 +36,6 @@ The corresponding call to warp is::
 
 The swirl transformation
 ````````````````````````
-
 Consider the coordinate :math:`(x, y)` in the output image.  The reverse
 mapping for the swirl transformation first computes, relative to a center
 :math:`(x_0, y_0)`, its polar coordinates,
@@ -60,9 +59,10 @@ and then transforms them according to
     \theta' = \phi + s \, e^{-\rho / r + \theta}
 
 where ``strength`` is a parameter for the amount of swirl, ``radius`` indicates
-the extent of the transform in pixels, and ``rotation`` adds a rotation angle.
-The transformation of ``radius`` into :math:`r` is to ensure that the
-transformation decays to :math:`\approx 1/1000^{\mathsf{th}}` within the specified radius.
+the swirl extent in pixels, and ``rotation`` adds a rotation angle.  The
+transformation of ``radius`` into :math:`r` is to ensure that the
+transformation decays to :math:`\approx 1/1000^{\mathsf{th}}` within the
+specified radius.
 """
 
 from skimage import data
diff --git a/skimage/transform/_swirl.py b/skimage/transform/_swirl.py
index 0e144ed1..6fd92c44 100644
--- a/skimage/transform/_swirl.py
+++ b/skimage/transform/_swirl.py
@@ -7,9 +7,12 @@ 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)
 
-    rho = np.sqrt((x - x0)**2 + (y - y0)**2)
     theta = rotation + strength * \
             np.exp(-rho / radius) + \
             np.arctan2(y - y0, x - x0)
@@ -32,8 +35,8 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
     strength : float
         The amount of swirling applied.
     radius : float
-        The extent of the swirling in pixels.  The effect dies out
-        rapidly beyond radius.
+        The extent of the swirl in pixels.  The effect dies out
+        rapidly beyond `radius`.
     rotation : float
         Additional rotation applied to the image.
 
@@ -47,10 +50,11 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
     output_shape : tuple or ndarray
         Size of the generated output image.
     order : int
-        Order of splines used in interpolation, passed as-is to ndimage.
+        Order of splines used in interpolation.  See
+        `scipy.ndimage.map_coordinates` for detail.
     mode : string
-        How to handle values outside the image borders, passed as-is
-        to ndimage.
+        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.
@@ -65,6 +69,6 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
                  'strength': strength,
                  'radius': radius}
 
-    return warp(image, _swirl_mapping, tf_args=warp_args,
+    return warp(image, _swirl_mapping, map_args=warp_args,
                 output_shape=output_shape,
                 order=order, mode=mode, cval=cval)
diff --git a/skimage/transform/_warp.py b/skimage/transform/_warp.py
index 00541fef..76f4fb5d 100644
--- a/skimage/transform/_warp.py
+++ b/skimage/transform/_warp.py
@@ -7,13 +7,21 @@ from skimage.util import img_as_float
 eps = np.finfo(float).eps
 
 def _stackcopy(a, b):
-    """a[:,:,0] = a[:,:,1] = ... = b"""
+    """Copy b into each color layer of a, such that::
+
+      a[:,:,0] = a[:,:,1] = ... = b
+
+    Notes
+    -----
+    Color images are stored as an ``MxNx3`` or ``MxNx4`` arrays.
+
+    """
     if a.ndim == 3:
         a.transpose().swapaxes(1, 2)[:] = b
     else:
         a[:] = b
 
-def warp(image, coord_tf, tf_args={},
+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.
 
@@ -21,20 +29,20 @@ def warp(image, coord_tf, tf_args={},
     ----------
     image : 2-D array
         Input image.
-    coord_tf : callable xy = f(xy, **kwargs)
-        Function that transforms an Nx2 array of ``(x, y)`` coordinates
-        in the *output image* into their corresponding coordinates in the
-        *source image*.  Note that this is a reverse mapping (also
-        see examples below).
-    tf_args : dict, optional
-        Keyword arguments passed to `coord_tf`.
+    reverse_map : callable xy = f(xy, **kwargs)
+        Reverse coordinate map.  A function that transforms an Nx2 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.
+        Order of splines used in interpolation. See
+        `scipy.ndimage.map_coordinates` for detail.
     mode : string
-        How to handle values outside the image borders.  Passed as-is
-        to ndimage.
+        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.
@@ -65,17 +73,24 @@ def warp(image, coord_tf, tf_args={},
 
     coords = np.empty(np.r_[3, output_shape], dtype=float)
 
-    # Construct transformed coordinates
+    ## 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
 
-    tf_coords = coord_tf(tf_coords, **tf_args)
+    # 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)
 
-    # y-coordinate mapping
+    # Place the y-coordinate mapping
     _stackcopy(coords[1, ...], tf_coords[0, ...])
 
-    # x-coordinate mapping
+    # Place the x-coordinate mapping
     _stackcopy(coords[0, ...], tf_coords[1, ...])
 
     # colour-coordinate mapping
diff --git a/skimage/transform/tests/test_swirl.py b/skimage/transform/tests/test_swirl.py
index e3fcc02e..d71f8231 100644
--- a/skimage/transform/tests/test_swirl.py
+++ b/skimage/transform/tests/test_swirl.py
@@ -6,11 +6,10 @@ 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'}
-    unswirled = tf.swirl(
-        tf.swirl(image, strength=10, **swirl_params),
-        strength=-10, **swirl_params
-        )
+    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