From 927157e7e508efbb0628f3e1f77bfb2bf60359ca Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 4 Oct 2009 01:34:30 +0200
Subject: [PATCH] Add projective transformation / homography.

---
 scikits/image/transform/__init__.py           |  2 +
 scikits/image/transform/project.py            | 99 +++++++++++++++++++
 scikits/image/transform/tests/test_project.py | 21 ++++
 3 files changed, 122 insertions(+)
 create mode 100644 scikits/image/transform/project.py
 create mode 100644 scikits/image/transform/tests/test_project.py

diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py
index 6eaf7dee..91711c6e 100644
--- a/scikits/image/transform/__init__.py
+++ b/scikits/image/transform/__init__.py
@@ -1,2 +1,4 @@
 from hough_transform import *
 from finite_radon_transform import *
+from project import *
+
diff --git a/scikits/image/transform/project.py b/scikits/image/transform/project.py
new file mode 100644
index 00000000..a6bc05df
--- /dev/null
+++ b/scikits/image/transform/project.py
@@ -0,0 +1,99 @@
+"""Image projection.
+
+"""
+
+import numpy as np
+from scipy.ndimage import interpolation as ndii
+
+__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.
+
+    For each pixel, given its homogeneous coordinate :math:`\mathbf{x}
+    = [x, y, 1]^T`, its target position is calculated by multiplying
+    with the given matrix, :math:`H`, to give :math:`H \mathbf{x}`.
+    E.g., to rotate by theta degrees clockwise, the matrix should be
+
+    [[cos(theta) -sin(theta) 0]
+     [sin(theta)  cos(theta) 0]
+     [0            0         1]]
+
+    or, to translate x by 10 and y by 20,
+
+    [[1 0 10]
+     [0 1 20]
+     [0 0 1 ]].
+
+    Parameters
+    ----------
+    image : 2-D array
+        Input image.
+    H : array of shape ``(3, 3)``
+        Transformation matrix H that defines the homography.
+    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.
+
+    """
+    if image.ndim < 2:
+        raise ValueError("Input must have more than 1 dimension.")
+
+    image = np.atleast_3d(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]
+    rows, cols = np.mgrid[:rows, :cols]
+    tf_coords = np.empty(shape=cols.shape,
+                         dtype=[('cols', float),
+                                ('rows', float),
+                                ('z', float)])
+    tf_coords['cols'], tf_coords['rows'] = cols, rows
+    tf_coords['z'] = 1
+    tf_coords = tf_coords.view((float, 3))
+
+    tf_coords = np.dot(tf_coords, np.linalg.inv(H).transpose())
+    tf_coords[np.absolute(tf_coords) < eps] = 0.
+
+    # normalize coordinates
+    tf_coords[..., :2] /= tf_coords[..., 2, np.newaxis]
+
+    # y-coordinate mapping
+    _stackcopy(coords[0,...], tf_coords[...,1])
+
+    # x-coordinate mapping
+    _stackcopy(coords[1,...], tf_coords[...,0])
+
+    # colour-coordinate mapping
+    coords[2,...] = range(bands)
+
+    # Prefilter not necessary for order 1 interpolation
+    prefilter = order > 1
+    mapped = ndii.map_coordinates(image, coords, prefilter=prefilter,
+                                  mode=mode, order=order, cval=cval)
+
+    return mapped.squeeze()
diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py
new file mode 100644
index 00000000..44c0e557
--- /dev/null
+++ b/scikits/image/transform/tests/test_project.py
@@ -0,0 +1,21 @@
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+
+from scikits.image.transform.project import _stackcopy, homography
+
+def test_stackcopy():
+    layers = 4
+    x = np.empty((3, 3, layers))
+    y = np.eye(3, 3)
+    _stackcopy(x, y)
+    for i in range(layers):
+        assert_array_almost_equal(x[...,i], y)
+
+def test_homography():
+    x = np.arange(9).reshape((3, 3)) + 1
+    theta = -np.pi/2
+    M = np.array([[np.cos(theta),-np.sin(theta),0],
+                  [np.sin(theta), np.cos(theta),2],
+                  [0,             0,            1]])
+    x90 = homography(x, M, order=1)
+    assert_array_almost_equal(x90, np.rot90(x))