ImProc ass2: implemented perspective transform.

improc/ass2/affine.py

 from interpolation import pv
-from pylab import array, matrix, zeros, show
-from numpy.linalg import lstsq, inv
+from pylab import array, matrix, zeros, show, dot, lstsq, inv, svd
 
 # url: http://www.leptonica.com/affine.html
 
@@ -35,26 +34,59 @@ def affine_transform(image, p1, p2, p3, width, height):
     # Construct the transformed image
     result = zeros((width, height))
 
-    for y in xrange(height):
-        for x in xrange(width):
+    for x in xrange(width):
+        for y in xrange(height):
             orig_pos = (A * array([[x, y, 1]]).T).T.tolist()[0]
             result[x][y] = pv(image, orig_pos[0], orig_pos[1], 'linear')
 
     return result
 
+def perspective_transform(image, p1, p2, p3, p4, width, height):
+    """Warps the parallelogram defined by four corner points onto a new image
+    of size MxN."""
+
+    x1, y1 = p1
+    x2, y2 = p2
+    x3, y3 = p3
+    x4, y4 = p4
+
+    M = array([[x1, y1, 1, 0, 0, 0,            0,            0,       0],
+               [0, 0, 0, x1, y1, 1,            0,            0,       0],
+               [x2, y2, 1, 0, 0, 0,  -width * x2,  -width * y2,  -width],
+               [0, 0, 0, x2, y2, 1,            0,            0,       0],
+               [x3, y3, 1, 0, 0, 0,            0,            0,       0],
+               [0, 0, 0, x3, y3, 1, -height * x3, -height * y3, -height],
+               [x4, y4, 1, 0, 0, 0,  -width * x4,  -width * y4,  -width],
+               [0, 0, 0, x4, y4, 1, -height * x4, -height * y4, -height]])
+
+    # Calculate transform matrix
+    A = inv(svd(M)[2].T[:,-1].reshape(3, 3))
+
+    # Construct the transformed image
+    result = zeros((width, height, 4))
+
+    for x in xrange(width):
+        for y in xrange(height):
+            p = dot(A, array([[x], [y], [1]]))
+            result[x][y] = pv(image, p[0] / p[2], p[1] / p[2], 'linear')
+
+    return result
+
 if __name__ == '__main__':
     from pylab import imread, imshow
 
-    image = imread('cameraman.png')
+    #image = imread('cameraman.png')
+    image = imread('flyeronground.png')
 
     half_width = image.shape[0] / 2
     half_height = image.shape[1] / 2
 
-    M, N = image.shape
+    M, N = image.shape[:2]
 
-    # TODO: validate MxN -> width x height
-    result = affine_transform(image, (0., half_height), (half_width, 0.),
-                              (image.shape[0], half_height), 128, 128)
+    #result = affine_transform(image, (0., half_height), (half_width, 0.),
+    #                          (image.shape[0], half_height), 128, 128)
+    result = perspective_transform(image, (149, 590), (101, 376), (301, 209),
+                                   (393, 392), 64, 64)
 
     imshow(result, cmap='gray')
     show()

improc/ass2/interpolation.py

@@ -12,10 +12,8 @@ def pv(image, x, y, method):
         return image[round(x)][round(y)] if in_image(image, x, y) else bgcolor
 
     if method == 'linear':
-        x1 = floor(x)
-        x2 = floor(x + 1)
-        y1 = floor(y)
-        y2 = floor(y + 1)
+        x1, y1 = floor(x), floor(y)
+        x2, y2 = floor(x + 1), floor(y + 1)
 
         if not in_image(image, x1, y1) or not in_image(image, x2, y2):
             return bgcolor