From c22a176e8d8c368e108cef65144e062abbd87121 Mon Sep 17 00:00:00 2001
From: James Bergstra <james.bergstra@gmail.com>
Date: Mon, 20 Aug 2012 19:13:04 -0400
Subject: [PATCH] ENH: moved code generating coords out of warp() to
 _warp_coords

---
 skimage/transform/_geometric.py | 63 ++++++++++++++++++++++-----------
 1 file changed, 42 insertions(+), 21 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 5baaeba1..9d394372 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -27,6 +27,45 @@ def _stackcopy(a, b):
         a[:] = b
 
 
+def _build_coords(orows, ocols, bands, coord_transform_fn,
+        dtype='float64'):
+    """
+    Return coords object suitable for scipy.ndimage.map_coordinates
+
+    Parameters
+    ----------
+    orows: number of output rows
+    ocols: number of output columns
+    bands: number of color bands (aka channels)
+    coord_transform_fn: something like GeometricTransform.inverse_map
+    interp_dtype: precision of interpolation e.g. 'float32' or 'float64'
+
+    """
+
+    coords = np.empty(np.r_[3, (orows, ocols, bands)], dtype=dtype)
+
+    # Reshape grid coordinates into a (P, 2) array of (x, y) pairs
+    tf_coords = np.indices((ocols, orows), dtype=dtype).reshape(2, -1).T
+
+    # Map each (x, y) pair to the source image according to
+    # the user-provided mapping
+    tf_coords = coord_transform_fn(tf_coords)
+
+    # Reshape back to a (2, M, N) coordinate grid
+    tf_coords = tf_coords.T.reshape((-1, ocols, orows)).swapaxes(1, 2)
+
+    # Place the y-coordinate mapping
+    _stackcopy(coords[1, ...], tf_coords[0, ...])
+
+    # Place the x-coordinate mapping
+    _stackcopy(coords[0, ...], tf_coords[1, ...])
+
+    # colour-coordinate mapping
+    coords[2, ...] = range(bands)
+
+    return coords
+
+
 class GeometricTransform(object):
     """Perform geometric transformations on a set of coordinates.
 
@@ -753,30 +792,12 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     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
+    def coord_transform_fn(*args):
+        return inverse_map(*args, **map_args)
 
-    # 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)
+    coords = _build_coords(rows, cols, bands, coord_transform_fn)
 
     # Prefilter not necessary for order 1 interpolation
     prefilter = order > 1