ENH: Implement fast coordinate transformations.

doc/examples/plot_swirl.py

@@ -0,0 +1,83 @@
+r"""
+=====
+Swirl
+=====
+
+Image swirling is a non-linear image deformation that creates a whirlpool
+effect.  This example describes the implementation of this transform in
+``skimage``, as well as the underlying warp mechanism.
+
+Image warping
+`````````````
+When applying a geometric transformation on an image, we typically make use of
+a reverse mapping, i.e., for each pixel in the output image, we compute its
+corresponding position in the input.  The reason is that, if we were to do it
+the other way around (map each input pixel to its new output position), some
+pixels in the output may be left empty.  On the other hand, each output
+coordinate has exactly one corresponding location in (or outside) the input
+image, and even if that position is non-integer, we may use interpolation to
+compute the corresponding image value.
+
+Performing a reverse mapping
+````````````````````````````
+To perform a geometric warp in ``skimage``, you simply need to provide the
+reverse mapping to the ``skimage.transform.warp`` function.  E.g., consider the
+case where we would like to shift an image 50 pixels to the left.  The reverse
+mapping for such a shift would be::
+
+    def shift_left(xy):
+        xy[:, 0] += 50
+        return xy
+
+The corresponding call to warp is::
+
+    from skimage.transform import warp
+    warp(image, shift_left)
+
+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,
+
+.. math::
+
+    \theta = \arctan(y/x)
+
+    \rho = \sqrt{(x - x_0)^2 + (y - y_0)^2},
+
+and then transforms them according to
+
+.. math::
+
+    r = \ln(2) \, \mathtt{radius} / 5
+
+    \phi = \mathtt{rotation}
+
+    s = \mathtt{strength}
+
+    \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.
+"""
+
+from skimage import data
+from skimage.transform import swirl
+
+import matplotlib.pyplot as plt
+
+image = data.checkerboard()
+swirled = swirl(image, rotation=0, strength=10, radius=120, order=2)
+
+f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3))
+
+ax0.imshow(image, cmap=plt.cm.gray, interpolation='none')
+ax0.axis('off')
+ax1.imshow(swirled, cmap=plt.cm.gray, interpolation='none')
+ax1.axis('off')
+
+plt.show()

skimage/transform/__init__.py

@@ -4,3 +4,5 @@ from .finite_radon_transform import *
 from .project import *
 from ._project import homography as fast_homography
 from .integral import *
+from ._warp import warp
+from ._swirl import swirl

skimage/transform/_swirl.py

@@ -0,0 +1,70 @@
+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
+    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)
+
+    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 swirling 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, passed as-is to ndimage.
+    mode : string
+        How to handle values outside the image borders, passed as-is
+        to ndimage.
+    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, tf_args=warp_args,
+                output_shape=output_shape,
+                order=order, mode=mode, cval=cval)

skimage/transform/_warp.py

@@ -0,0 +1,91 @@
+__all__ = ['warp']
+
+import numpy as np
+from scipy import ndimage
+from skimage.util import img_as_float
+
+eps = np.finfo(float).eps
+
+def _stackcopy(a, b):
+    """a[:,:,0] = a[:,:,1] = ... = b"""
+    if a.ndim == 3:
+        a.transpose().swapaxes(1, 2)[:] = b
+    else:
+        a[:] = b
+
+def warp(image, coord_tf, tf_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.
+    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`.
+    output_shape : tuple (rows, cols)
+        Shape of the output image generated.
+    order : int
+        Order of splines used in interpolation.
+    mode : string
+        How to handle values outside the image borders.  Passed as-is
+        to ndimage.
+    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]
+    tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T
+
+    tf_coords = coord_tf(tf_coords, **tf_args)
+    tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
+
+    # y-coordinate mapping
+    _stackcopy(coords[1, ...], tf_coords[0, ...])
+
+    # 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)

skimage/transform/project.py

@@ -4,18 +4,12 @@
 
 import numpy as np
 from scipy.ndimage import interpolation as ndii
+from .warp import _stackcopy
 
 __all__ = ['homography']
 
 eps = np.finfo(float).eps
 
-def _stackcopy(a, b):
-    """a[:,:,0] = a[:,:,1] = ... = b"""
-    if a.ndim == 3:
-        a.transpose().swapaxes(1, 2)[:] = b
-    else:
-        a[:] = b
-
 def homography(image, H, output_shape=None, order=1,
                mode='constant', cval=0.):
     """Perform a projective transformation (homography) on an image.
@@ -106,6 +100,8 @@ def homography(image, H, output_shape=None, order=1,
 
     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]

skimage/transform/tests/test_project.py

@@ -1,7 +1,7 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal
 
-from skimage.transform.project import _stackcopy
+from skimage.transform._warp import _stackcopy
 from skimage.transform import homography, fast_homography
 from skimage import data
 from skimage.color import rgb2gray

skimage/transform/tests/test_swirl.py

@@ -0,0 +1,18 @@
+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'}
+    unswirled = tf.swirl(
+        tf.swirl(image, strength=10, **swirl_params),
+        strength=-10, **swirl_params
+        )
+
+    assert np.mean(np.abs(image - unswirled)) < 0.01
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()