speedup from_film() projection.

1 files changed, 23 insertions, 14 deletions

diff --git a/project.py b/project.py
index 0e3626f..10254f2 100644
--- a/project.py
+++ b/project.py
@@ -1,28 +1,37 @@
 import numpy as np
 from itertools import repeat
-from chroma.transform import rotate
+from chroma.transform import rotate, normalize
+
+def from_film(position, axis1=(0,0,1), axis2=(1,0,0), size=(800,600),
+              width=0.035, focal_length=0.018):
+    """Project rays from a piece of film located at whose focal point is
+    located at `position`. `axis1` and `axis2` specify the vectors pointing
+    along the length and height of the film respectively.
+    """
 
-def from_film(position, axis1=(0,0,1), axis2=(1,0,0), size=(800,600), width=0.035, focal_length=0.018):
     height = width*(size[1]/float(size[0]))
 
-    x = np.linspace(width/2, -width/2, size[0])
-    z = np.linspace(-height/2, height/2, size[1])
+    axis1 = normalize(axis1)
+    axis2 = normalize(axis2)
 
-    zz, xx = np.meshgrid(z,x)
+    dx0 = width/size[0]
+    dx1 = height/size[1]
 
-    grid = np.array(zip(xx.flatten(),np.zeros(xx.size),zz.flatten()))
+    # for i in range(size[0]):
+    #     for j in range(size[1]):
+    #         grid[i*size[1]+j] = axis1*dx1*j - axis2*dx0*i
 
-    axis1 = np.asarray(axis1)/np.linalg.norm(axis1)
-    axis2 = np.asarray(axis2)/np.linalg.norm(axis2)
+    x = np.arange(size[0])
+    y = np.arange(size[1])
 
-    assert axis1.dot(axis2) == 0.0
+    yy, xx = np.meshgrid(y,x)
 
-    if (axis1 != (0,0,1)).any():
-        grid = rotate(grid, np.arccos(axis1.dot((0,0,1))), np.cross(axis1,(0,0,1)))
+    n = size[0]*size[1]
 
-    if (axis2 != (1,0,0)).any():
-        grid = rotate(grid, np.arccos(axis2.dot((1,0,0))), np.cross(axis2,(1,0,0)))
+    grid = -np.tile(axis2, (n,1))*xx.ravel()[:,np.newaxis]*dx0 + \
+        np.tile(axis1, (n,1))*yy.ravel()[:,np.newaxis]*dx1
 
+    grid += axis2*width/2 - axis1*height/2
     grid -= np.cross(axis1,axis2)*focal_length
 
-    return grid+position, -grid/np.apply_along_axis(np.linalg.norm,1,grid)[:,np.newaxis]
+    return grid+position, normalize(-grid)