From 2d94a273f6c77aa0c2d0f70092fd8966f141b068 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 17:22:34 -0700
Subject: [PATCH] Fix mirror mode. Update docstrings. Update tests

---
 scikits/image/transform/_project.pyx          | 61 ++++++++++++++++---
 scikits/image/transform/tests/test_project.py | 36 +++++++++++
 2 files changed, 88 insertions(+), 9 deletions(-)

diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx
index e2d08a02..d1788ab8 100644
--- a/scikits/image/transform/_project.pyx
+++ b/scikits/image/transform/_project.pyx
@@ -1,5 +1,7 @@
 #cython: cdivison=True boundscheck=False
 
+__all__ = ['homography']
+
 cimport cython
 cimport numpy as np
 
@@ -16,6 +18,22 @@ cdef extern from "math.h":
 
 cdef double get_pixel(double *image, int rows, int cols,
                       int r, int c, char mode, double cval=0):
+    """Get a pixel from the image, taking wrapping mode into consideration.
+
+    Parameters
+    ----------
+    image : *double
+        Input image.
+    rows, cols : int
+        Dimensions of image.
+    r, c : int
+        Position at which to get the pixel.
+    mode : {'C', 'W', 'M'}
+        Wrapping mode.  Constant, Wrap or Mirror.
+    cval : double
+        Constant value to use for mode constant.
+    
+    """
     if mode == 'C':
         if (r < 0) or (r > cols - 1) or (c < 0) or (c > cols - 1):
             return cval
@@ -26,10 +44,28 @@ cdef double get_pixel(double *image, int rows, int cols,
                      coord_map(cols, c, mode)]
 
 cdef int coord_map(int dim, int coord, char mode):
+    """
+    Wrap a coordinate, according to a given dimension and mode.
+    
+    Parameters
+    ----------
+    dim : int
+        Maximum coordinate.
+    coord : int
+        Coord provided by user.  May be < 0 or > dim.
+    mode : {'W', 'M'}
+        Whether to wrap or mirror the coordinate if it
+        falls outside [0, dim).
+    
+    """
     dim = dim - 1
     if mode == 'M': # mirror
         if (coord < 0):
-            return <int>(-coord % dim)
+            # How many times times does the coordinate wrap?
+            if (<int>(-coord / dim) % 2 != 0):
+                return dim - <int>(-coord % dim)
+            else:
+                return <int>(-coord % dim)
         elif (coord > dim):
             if (<int>(coord / dim) % 2 != 0):
                 return <int>(dim - (coord % dim))
@@ -44,13 +80,19 @@ cdef int coord_map(int dim, int coord, char mode):
     return coord
 
 cdef tf(double x, double y, double* H, double *x_, double *y_):
-    cdef double xx, yy, zz
+    """Apply a homography to a coordinate.
 
-    ## print
-    ## print H[0], H[1], H[2]
-    ## print H[3], H[4], H[5]
-    ## print H[6], H[7], H[8]
-    ## print
+    Parameters
+    ----------
+    x, y : double
+        Input coordinate.
+    H : (3,3) *double
+        Transformation matrix.
+    x_, y_ : *double
+        Output coordinate.
+
+    """
+    cdef double xx, yy, zz
 
     xx = H[0] * x + H[1] * y + H[2]
     yy = H[3] * x + H[4] * y + H[5]
@@ -108,7 +150,8 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
 
     """
 
-    cdef np.ndarray[dtype=np.double_t, ndim=2] img = image
+    cdef np.ndarray[dtype=np.double_t, ndim=2] img = \
+         np.asarray(image, dtype=np.double)
     cdef np.ndarray[dtype=np.double_t, ndim=2] M = \
          np.ascontiguousarray(np.linalg.inv(H))
 
@@ -134,7 +177,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
     columns = img.shape[1]
 
     cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
-         np.zeros((out_r, out_c), dtype=np.float64)
+         np.zeros((out_r, out_c), dtype=np.double)
     
     cdef int tfr, tfc, r_int, c_int
     cdef double y0, y1, y2, y3
diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py
index 44c0e557..7ae24049 100644
--- a/scikits/image/transform/tests/test_project.py
+++ b/scikits/image/transform/tests/test_project.py
@@ -2,6 +2,8 @@ import numpy as np
 from numpy.testing import assert_array_almost_equal
 
 from scikits.image.transform.project import _stackcopy, homography
+from scikits.image.transform._project import homography as fast_homography
+from scikits.image import data
 
 def test_stackcopy():
     layers = 4
@@ -19,3 +21,37 @@ def test_homography():
                   [0,             0,            1]])
     x90 = homography(x, M, order=1)
     assert_array_almost_equal(x90, np.rot90(x))
+
+def test_fast_homography():
+    img = data.lena()
+    img = img[:, :100]
+    
+    theta = np.deg2rad(30)
+    scale = 0.5
+    tx, ty = 50, 50
+
+    H = np.eye(3)
+    S = scale * np.sin(theta)
+    C = scale * np.cos(theta)
+
+    H[:2, :2] = [[C, -S], [S, C]]
+    H[:2, 2] = [tx, ty]
+
+    for mode in ('constant', 'mirror', 'wrap'):
+        print 'Transform mode:', mode
+
+        p0 = homography(img, H, mode=mode)
+        p1 = fast_homography(img, H, mode=mode)
+        p1 = np.round(p1)
+
+        ## import matplotlib.pyplot as plt
+        ## plt.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
+        ## plt.show()
+
+        d = np.mean(np.abs(p0 - p1))
+        assert d < 0.1
+    
+
+if __name__ == "__main__":
+    from numpy.testing import run_module_suite
+    run_module_suite()