Merge pull request #253 from jaberg/warp_cval_outside_unit_interval

skimage/transform/__init__.py

@@ -3,7 +3,7 @@ from .radon_transform import *
 from .finite_radon_transform import *
 from ._project import homography as fast_homography
 from .integral import *
-from ._geometric import (warp, estimate_transform,
+from ._geometric import (warp, warp_coords, estimate_transform,
                          SimilarityTransform, AffineTransform,
                          ProjectiveTransform, PolynomialTransform)
 from ._warps import swirl, homography

skimage/transform/_geometric.py

@@ -698,6 +698,83 @@ def matrix_transform(coords, matrix):
     return ProjectiveTransform(matrix)(coords)
 
 
+def warp_coords(orows, ocols, bands, coord_transform_fn,
+                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
+        Return input coordinates for given output coordinates
+    dtype : np.dtype or string
+        dtype for return value (sane choices: float32 or float64)
+
+    Returns
+    -------
+    coords : (3, orows, ocols, 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`.
+
+    Notes
+    -----
+    This is a lower-level routine that produces the source coordinates used by
+    `warp()`.
+
+    It is provided separately from `warp` to give additional flexibility to
+    users who would like, for example, to re-use a particular coordinate
+    mapping, to use specific dtypes at various points along the the
+    image-warping process, or to implement different post-processing logic
+    than `warp` performs after the call to `ndimage.map_coordinates`.
+
+
+    Examples
+    --------
+    Produce a coordinate map that Shifts an image to the right:
+
+    >>> from skimage import data
+    >>> from scipy.ndimage import map_coordinates
+    >>>
+    >>> def shift_right(xy):
+    ...     xy[:, 0] -= 10
+    ...     return xy
+    >>>
+    >>> coords = warp_coords(30, 30, 3, shift_right)
+    >>> image = data.lena().astype(np.float32)
+    >>> warped_image = map_coordinates(image, coords)
+
+    """
+
+    coords = np.empty((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
+
+
 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.
@@ -721,7 +798,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     mode : string
         How to handle values outside the image borders.  See
         `scipy.ndimage.map_coordinates` for detail.
-    cval : string
+    cval : float
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
 
@@ -753,30 +830,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 = warp_coords(rows, cols, bands, coord_transform_fn)
 
     # Prefilter not necessary for order 1 interpolation
     prefilter = order > 1
@@ -785,4 +844,10 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
     # The spline filters sometimes return results outside [0, 1],
     # so clip to ensure valid data
-    return np.clip(mapped.squeeze(), 0, 1)
+    clipped = np.clip(mapped, 0, 1)
+
+    if mode == 'constant' and not (0 <= cval <= 1):
+        clipped[mapped == cval] = cval
+
+    # Remove singleton dim introduced by atleast_3d
+    return clipped.squeeze()

skimage/transform/tests/test_warps.py

@@ -1,8 +1,11 @@
 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, homography, fast_homography,
-                               SimilarityTransform, ProjectiveTransform)
+from skimage.transform import (warp, warp_coords, fast_homography,
+                               AffineTransform,
+                               ProjectiveTransform,
+                               SimilarityTransform)
 from skimage import transform as tf, data, img_as_float
 from skimage.color import rgb2gray
 
@@ -25,16 +28,19 @@ 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)
+    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 = warp(x,
+            inverse_map=ProjectiveTransform(M).inverse,
+            order=1)
     assert_array_almost_equal(x90, np.rot90(x))
 
 
 def test_fast_homography():
-    img = rgb2gray(data.lena())
+    img = rgb2gray(data.lena()).astype(np.uint8)
     img = img[:, :100]
 
     theta = np.deg2rad(30)
@@ -74,5 +80,33 @@ def test_swirl():
     assert np.mean(np.abs(image - unswirled)) < 0.01
 
 
+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)
+
+
+def test_warp_identity():
+    lena = img_as_float(rgb2gray(data.lena()))
+    assert len(lena.shape) == 2
+    assert np.allclose(lena, warp(lena, AffineTransform(rotation=0)))
+    assert not np.allclose(lena, warp(lena, AffineTransform(rotation=0.1)))
+    rgb_lena = np.transpose(np.asarray([lena, np.zeros_like(lena), lena]),
+                            (1, 2, 0))
+    warped_rgb_lena = warp(rgb_lena, AffineTransform(rotation=0.1))
+    assert np.allclose(rgb_lena, warp(rgb_lena, AffineTransform(rotation=0)))
+    assert not np.allclose(rgb_lena, warped_rgb_lena)
+    # assert no cross-talk between bands
+    assert np.all(0 == warped_rgb_lena[:, :, 1])
+
+
+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)
+    warped_image1 = map_coordinates(image[:, :, 0], coords[:2])
+
+
 if __name__ == "__main__":
     run_module_suite()