From 16736eeae08d9627ad751c65919900ae9191a08f Mon Sep 17 00:00:00 2001
From: Anthony LaTorre <telatorre@gmail.com>
Date: Sat, 7 May 2011 20:32:35 -0400
Subject: tie fighter

---
 intersect.cu                |  72 --------------
 linalg.h                    | 116 -----------------------
 matrix.h                    | 224 -------------------------------------------
 models/40mmcube.stl         | Bin 0 -> 684 bytes
 models/ccube4.STL           | Bin 0 -> 472984 bytes
 models/companioncube.stl    | Bin 0 -> 862334 bytes
 models/lionsolid.stl        | Bin 0 -> 3717984 bytes
 models/tie_interceptor6.stl | Bin 0 -> 442284 bytes
 rotate.h                    |  24 -----
 src/intersect.cu            | 118 +++++++++++++++++++++++
 src/linalg.h                | 116 +++++++++++++++++++++++
 src/matrix.h                | 225 ++++++++++++++++++++++++++++++++++++++++++++
 src/rotate.h                |  24 +++++
 stl.py                      |  25 ++++-
 test.py                     |  88 +++++++++++++++++
 15 files changed, 591 insertions(+), 441 deletions(-)
 delete mode 100644 intersect.cu
 delete mode 100644 linalg.h
 delete mode 100644 matrix.h
 create mode 100644 models/40mmcube.stl
 create mode 100644 models/ccube4.STL
 create mode 100644 models/companioncube.stl
 create mode 100644 models/lionsolid.stl
 create mode 100644 models/tie_interceptor6.stl
 delete mode 100644 rotate.h
 create mode 100644 src/intersect.cu
 create mode 100644 src/linalg.h
 create mode 100644 src/matrix.h
 create mode 100644 src/rotate.h
 create mode 100644 test.py

diff --git a/intersect.cu b/intersect.cu
deleted file mode 100644
index d7fdaa5..0000000
--- a/intersect.cu
+++ /dev/null
@@ -1,72 +0,0 @@
-//-*-c-*-
-
-__device__ __host__ Matrix inv(const Matrix&m, float& determinant)
-{
-	determinant = det(m);
-
-	return make_matrix(m.a11*m.a22 - m.a12*m.a21,
-			   m.a02*m.a21 - m.a01*m.a22,
-			   m.a01*m.a12 - m.a02*m.a11,
-			   m.a12*m.a20 - m.a10*m.a22,
-			   m.a00*m.a22 - m.a02*m.a20,
-			   m.a02*m.a10 - m.a00*m.a12,
-			   m.a10*m.a21 - m.a11*m.a20,
-			   m.a01*m.a20 - m.a00*m.a21,
-			   m.a00*m.a11 - m.a01*m.a10)/determinant;
-}
-
-__device__ bool intersect_triangle(const float3 &x, const float3 &p, float3 *vertex, float3 &intersection)
-{
-	float determinant;
-
-	float3 u = inv(make_matrix(vertex[0]-vertex[2],vertex[1]-vertex[2],-p), determinant)*(x-vertex[2]);
-
-	if (determinant == 0.0)
-		return false;
-
-	if (u.x < 0.0 || u.y < 0.0 || u.z < 0.0 || (1-u.x-u.y) < 0.0)
-		return false;
-
-	intersection = x + p*u.z;
-
-	return true;
-}
-
-extern "C"
-{
-
-__global__ void intersect_triangle_mesh(float3 *x, float3 *p, int n, float3 *mesh, int *state, float3 *w)
-{
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
-
-	bool hit = false;
-
-	float distance, min_distance;
-	float3 intersection, min_intersection;
-
-	int i;
-	for (i=0; i < n; i++)
-	{
-		if (intersect_triangle(x[idx], p[idx], mesh+3*i, intersection))
-		{
-			if (!hit)
-			{
-				state[idx] = 1;
-				min_distance = norm(intersection-x[idx]);
-				min_intersection = intersection;
-				continue;
-			}
-
-			distance = norm(w[idx]-x[idx]);
-
-			if (distance < min_distance)
-			{
-				state[idx] = 1;
-				min_distance = distance;
-				min_intersection = intersection;
-			}
-		}
-	}
-}
-
-} // extern "c"
diff --git a/linalg.h b/linalg.h
deleted file mode 100644
index 593aa9d..0000000
--- a/linalg.h
+++ /dev/null
@@ -1,116 +0,0 @@
-#ifndef __LINALG_H__
-#define __LINALG_H__
-
-__device__ __host__ float3 operator- (const float3 &a)
-{
-	return make_float3(-a.x, -a.y, -a.z);
-}
-
-__device__ __host__ float3 operator+ (const float3 &a, const float3 &b)
-{
-	return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
-}
-
-__device__ __host__ void operator+= (float3 &a, const float3 &b)
-{
-	a.x += b.x;
-	a.y += b.y;
-	a.z += b.z;
-}
-
-__device__ __host__ float3 operator- (const float3 &a, const float3 &b)
-{
-	return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
-}
-
-__device__ __host__ void operator-= (float3 &a, const float3 &b)
-{
-	a.x -= b.x;
-	a.y -= b.y;
-	a.z -= b.z;
-}
-
-__device__ __host__ float3 operator+ (const float3 &a, const float &c)
-{
-	return make_float3(a.x+c, a.y+c, a.z+c);
-}
-
-__device__ __host__ void operator+= (float3 &a, const float &c)
-{
-	a.x += c;
-	a.y += c;
-	a.z += c;
-}
-
-__device__ __host__ float3 operator+ (const float &c, const float3 &a)
-{
-	return make_float3(c+a.x, c+a.y, c+a.z);
-}
-
-__device__ __host__ float3 operator- (const float3 &a, const float &c)
-{
-	return make_float3(a.x-c, a.y-c, a.z-c);
-}
-
-__device__ __host__ void operator-= (float3 &a, const float &c)
-{
-	a.x -= c;
-	a.y -= c;
-	a.z -= c;
-}
-
-__device__ __host__ float3 operator- (const float &c, const float3& a)
-{
-	return make_float3(c-a.x, c-a.y, c-a.z);
-}
-
-__device__ __host__ float3 operator* (const float3 &a, const float &c)
-{
-	return make_float3(a.x*c, a.y*c, a.z*c);
-}
-
-__device__ __host__ void operator*= (float3 &a, const float &c)
-{
-	a.x *= c;
-	a.y *= c;
-	a.z *= c;
-}
-
-__device__ __host__ float3 operator* (const float &c, const float3& a)
-{
-	return make_float3(c*a.x, c*a.y, c*a.z);
-}
-
-__device__ __host__ float3 operator/ (const float3 &a, const float &c)
-{
-	return make_float3(a.x/c, a.y/c, a.z/c);
-}
-
-__device__ __host__ void operator/= (float3 &a, const float &c)
-{
-	a.x /= c;
-	a.y /= c;
-	a.z /= c;
-}
-
-__device__ __host__ float3 operator/ (const float &c, const float3 &a)
-{
-	return make_float3(c/a.x, c/a.y, c/a.z);
-}
-
-__device__ __host__ float dot(const float3 &a, const float3 &b)
-{
-	return a.x*b.x + a.y*b.y + a.z*b.z;
-}
-
-__device__ __host__ float3 cross(const float3 &a, const float3 &b)
-{
-	return make_float3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
-}
-
-__device__ __host__ float norm(const float3 &a)
-{
-	return sqrtf(dot(a,a));
-}
-
-#endif
diff --git a/matrix.h b/matrix.h
deleted file mode 100644
index 1de8e8c..0000000
--- a/matrix.h
+++ /dev/null
@@ -1,224 +0,0 @@
-#ifndef __MATRIX_H__
-#define __MATRIX_H__
-
-struct Matrix
-{
-	float a00, a01, a02, a10, a11, a12, a20, a21, a22;
-};
-
-__device__ __host__ Matrix make_matrix(float a00, float a01, float a02,
-				       float a10, float a11, float a12,
-				       float a20, float a21, float a22)
-{
-	Matrix m = {a00, a01, a02, a10, a11, a12, a20, a21, a22};
-	return m;
-}
-
-__device__ __host__ Matrix make_matrix(float3 &u1, float3 &u2, float3 &u3)
-{
-	Matrix m = {u1.x, u2.x, u3.x, u1.y, u2.y, u3.y, u1.z, u2.z, u3.z};
-}
-
-__device__ __host__ Matrix operator- (const Matrix &m)
-{
-	return make_matrix(-m.a00, -m.a01, -m.a02,
-			   -m.a10, -m.a11, -m.a12,
-			   -m.a20, -m.a21, -m.a22);
-}
-
-__device__ __host__ float3 operator* (const Matrix &m, const float3 &a)
-{
-	return make_float3(m.a00*a.x + m.a01*a.y + m.a02*a.z,
-			   m.a10*a.x + m.a11*a.y + m.a12*a.z,
-			   m.a20*a.x + m.a21*a.y + m.a22*a.z);
-}
-
-__device__ __host__ Matrix operator+ (const Matrix &m, const Matrix &n)
-{
-	return make_matrix(m.a00+n.a00, m.a01+n.a01, m.a02+n.a02,
-			   m.a10+n.a10, m.a11+n.a11, m.a12+n.a12,
-			   m.a20+n.a20, m.a21+n.a21, m.a22+n.a22);
-}
-
-__device__ __host__ void operator+= (Matrix &m, const Matrix &n)
-{
-	m.a00 += n.a00;
-	m.a01 += n.a01;
-	m.a02 += n.a02;
-	m.a10 += n.a10;
-	m.a11 += n.a11;
-	m.a12 += n.a12;
-	m.a20 += n.a20;
-	m.a21 += n.a21;
-	m.a22 += n.a22;
-}
-
-__device__ __host__ Matrix operator- (const Matrix &m, const Matrix &n)
-{
-	return make_matrix(m.a00-n.a00, m.a01-n.a01, m.a02-n.a02,
-			   m.a10-n.a10, m.a11-n.a11, m.a12-n.a12,
-			   m.a20-n.a20, m.a21-n.a21, m.a22-n.a22);
-}
-
-__device__ __host__ void operator-= (Matrix &m, const Matrix& n)
-{
-	m.a00 -= n.a00;
-	m.a01 -= n.a01;
-	m.a02 -= n.a02;
-	m.a10 -= n.a10;
-	m.a11 -= n.a11;
-	m.a12 -= n.a12;
-	m.a20 -= n.a20;
-	m.a21 -= n.a21;
-	m.a22 -= n.a22;
-}
-
-__device__ __host__ Matrix operator* (const Matrix &m, const Matrix &n)
-{
-	return make_matrix(m.a00*n.a00 + m.a01*n.a10 + m.a02*n.a20,
-			   m.a00*n.a01 + m.a01*n.a11 + m.a02*n.a21,
-			   m.a00*n.a02 + m.a01*n.a12 + m.a02*n.a22,
-			   m.a10*n.a00 + m.a11*n.a10 + m.a12*n.a20,
-			   m.a10*n.a01 + m.a11*n.a11 + m.a12*n.a21,
-			   m.a10*n.a02 + m.a11*n.a12 + m.a12*n.a22,
-			   m.a20*n.a00 + m.a21*n.a10 + m.a22*n.a20,
-			   m.a20*n.a01 + m.a21*n.a11 + m.a22*n.a21,
-			   m.a20*n.a02 + m.a21*n.a12 + m.a22*n.a22);
-}
-
-__device__ __host__ Matrix operator+ (const Matrix &m, const float &c)
-{
-	return make_matrix(m.a00+c, m.a01+c, m.a02+c,
-			   m.a10+c, m.a11+c, m.a12+c,
-			   m.a20+c, m.a21+c, m.a22+c);
-}
-
-__device__ __host__ void operator+= (Matrix &m, const float &c)
-{
-	m.a00 += c;
-	m.a01 += c;
-	m.a02 += c;
-	m.a10 += c;
-	m.a11 += c;
-	m.a12 += c;
-	m.a20 += c;
-	m.a21 += c;
-	m.a22 += c;
-}
-
-__device__ __host__ Matrix operator+ (const float &c, const Matrix &m)
-{
-	return make_matrix(c+m.a00, c+m.a01, c+m.a02,
-			   c+m.a10, c+m.a11, c+m.a12,
-			   c+m.a20, c+m.a21, c+m.a22);
-}
-
-__device__ __host__ Matrix operator- (const Matrix &m, const float &c)
-{
-	return make_matrix(m.a00-c, m.a01-c, m.a02-c,
-			   m.a10-c, m.a11-c, m.a12-c,
-			   m.a20-c, m.a21-c, m.a22-c);
-}
-
-__device__ __host__ void operator-= (Matrix &m, const float &c)
-{
-	m.a00 -= c;
-	m.a01 -= c;
-	m.a02 -= c;
-	m.a10 -= c;
-	m.a11 -= c;
-	m.a12 -= c;
-	m.a20 -= c;
-	m.a21 -= c;
-	m.a22 -= c;
-}
-
-__device__ __host__ Matrix operator- (const float &c, const Matrix &m)
-{
-	return make_matrix(c-m.a00, c-m.a01, c-m.a02,
-			   c-m.a10, c-m.a11, c-m.a12,
-			   c-m.a20, c-m.a21, c-m.a22);
-}
-
-__device__ __host__ Matrix operator* (const Matrix &m, const float &c)
-{
-	return make_matrix(m.a00*c, m.a01*c, m.a02*c,
-			   m.a10*c, m.a11*c, m.a12*c,
-			   m.a20*c, m.a21*c, m.a22*c);
-}
-
-__device__ __host__ void operator*= (Matrix &m, const float &c)
-{
-	m.a00 *= c;
-	m.a01 *= c;
-	m.a02 *= c;
-	m.a10 *= c;
-	m.a11 *= c;
-	m.a12 *= c;
-	m.a20 *= c;
-	m.a21 *= c;
-	m.a22 *= c;
-}
-
-__device__ __host__ Matrix operator* (const float &c, const Matrix &m)
-{
-	return make_matrix(c*m.a00, c*m.a01, c*m.a02,
-			   c*m.a10, c*m.a11, c*m.a12,
-			   c*m.a20, c*m.a21, c*m.a22);
-}
-
-__device__ __host__ Matrix operator/ (const Matrix &m, const float &c)
-{
-	return make_matrix(m.a00/c, m.a01/c, m.a02/c,
-			   m.a10/c, m.a11/c, m.a12/c,
-			   m.a20/c, m.a21/c, m.a22/c);
-}
-
-__device__ __host__ void operator/= (Matrix &m, const float &c)
-{
-	m.a00 /= c;
-	m.a01 /= c;
-	m.a02 /= c;
-	m.a10 /= c;
-	m.a11 /= c;
-	m.a12 /= c;
-	m.a20 /= c;
-	m.a21 /= c;
-	m.a22 /= c;
-}
-
-__device__ __host__ Matrix operator/ (const float &c, const Matrix &m)
-{
-	return make_matrix(c/m.a00, c/m.a01, c/m.a02,
-			   c/m.a10, c/m.a11, c/m.a12,
-			   c/m.a20, c/m.a21, c/m.a22);
-}
-
-__device__ __host__ float det(const Matrix &m)
-{
-	return m.a00*(m.a11*m.a22 - m.a12*m.a21) - \
-		m.a10*(m.a01*m.a22 - m.a02*m.a21) + \
-		m.a20*(m.a01*m.a12 - m.a02*m.a11);
-}
-
-__device__ __host__ Matrix inv(const Matrix &m)
-{
-	return make_matrix(m.a11*m.a22 - m.a12*m.a21,
-			   m.a02*m.a21 - m.a01*m.a22,
-			   m.a01*m.a12 - m.a02*m.a11,
-			   m.a12*m.a20 - m.a10*m.a22,
-			   m.a00*m.a22 - m.a02*m.a20,
-			   m.a02*m.a10 - m.a00*m.a12,
-			   m.a10*m.a21 - m.a11*m.a20,
-			   m.a01*m.a20 - m.a00*m.a21,
-			   m.a00*m.a11 - m.a01*m.a10)/det(m);
-}
-
-__device__ __host__ Matrix outer(const float3 &a, const float3 &b)
-{
-	return make_matrix(a.x*b.x, a.x*b.y, a.x*b.z,
-			   a.y*b.x, a.y*b.y, a.y*b.z,
-			   a.z*b.x, a.z*b.y, a.z*b.z);
-}
-
-#endif
diff --git a/models/40mmcube.stl b/models/40mmcube.stl
new file mode 100644
index 0000000..779ff9c
Binary files /dev/null and b/models/40mmcube.stl differ
diff --git a/models/ccube4.STL b/models/ccube4.STL
new file mode 100644
index 0000000..126dc86
Binary files /dev/null and b/models/ccube4.STL differ
diff --git a/models/companioncube.stl b/models/companioncube.stl
new file mode 100644
index 0000000..09ea5a1
Binary files /dev/null and b/models/companioncube.stl differ
diff --git a/models/lionsolid.stl b/models/lionsolid.stl
new file mode 100644
index 0000000..e133e8f
Binary files /dev/null and b/models/lionsolid.stl differ
diff --git a/models/tie_interceptor6.stl b/models/tie_interceptor6.stl
new file mode 100644
index 0000000..970e8da
Binary files /dev/null and b/models/tie_interceptor6.stl differ
diff --git a/rotate.h b/rotate.h
deleted file mode 100644
index 52d6d6a..0000000
--- a/rotate.h
+++ /dev/null
@@ -1,24 +0,0 @@
-#ifndef __ROTATE_H__
-#define __ROTATE_H__
-
-__device__ const Matrix IDENTITY_MATRIX = {1,0,0,0,1,0,0,0,1};
-__device__ const Matrix ZERO_MATRIX = {0,0,0,0,0,0,0,0,0};
-
-__device__ __host__ Matrix make_rotation_matrix(float phi, const float3 &n)
-{
-	/* rotate points counterclockwise, when looking towards +infinity,
-	   through an angle `phi` about the axis `n`. */
-
-	float cos_phi = cosf(phi);
-	float sin_phi = sinf(phi);
-
-	return IDENTITY_MATRIX*cos_phi + (1-cos_phi)*outer(n,n) +
-		sin_phi*make_matrix(0,n.z,-n.y,-n.z,0,n.x,n.y,-n.x,0);
-}
-
-__device__ __host__ float3 rotate(const float3 &a, float phi, const float3 &n)
-{
-	return make_rotation_matrix(phi,n)*a;
-}
-
-#endif
diff --git a/src/intersect.cu b/src/intersect.cu
new file mode 100644
index 0000000..844f6fd
--- /dev/null
+++ b/src/intersect.cu
@@ -0,0 +1,118 @@
+//-*-c-*-
+
+__device__ __host__ Matrix inv(const Matrix&m, float& determinant)
+{
+	determinant = det(m);
+
+	return make_matrix(m.a11*m.a22 - m.a12*m.a21,
+			   m.a02*m.a21 - m.a01*m.a22,
+			   m.a01*m.a12 - m.a02*m.a11,
+			   m.a12*m.a20 - m.a10*m.a22,
+			   m.a00*m.a22 - m.a02*m.a20,
+			   m.a02*m.a10 - m.a00*m.a12,
+			   m.a10*m.a21 - m.a11*m.a20,
+			   m.a01*m.a20 - m.a00*m.a21,
+			   m.a00*m.a11 - m.a01*m.a10)/determinant;
+}
+
+__device__ bool intersect_triangle(const float3 &x, const float3 &p, float3 *vertex, float3 &intersection)
+{
+	float determinant;
+	float3 u = inv(make_matrix(vertex[1]-vertex[0],vertex[2]-vertex[0],-p), determinant)*(x-vertex[0]);
+
+	if (determinant == 0.0)
+		return false;
+
+	if (u.x < 0.0 || u.y < 0.0 || u.z < 0.0 || (1-u.x-u.y) < 0.0)
+		return false;
+
+	intersection = x + p*u.z;
+
+	return true;
+}
+
+__device__ int get_color(const float3 &p, float3 *vertex)
+{
+	float3 v1 = vertex[1] - vertex[0];
+	float3 v2 = vertex[2] - vertex[0];
+
+	float3 normal = cross(v1,v2);
+
+	float scale;
+	scale = dot(normal,-p)/(norm(normal)*norm(p));
+	if (scale < 0.0)
+		scale = dot(-normal,-p)/(norm(normal)*norm(p));
+
+	int rgb = floorf(255*scale);
+
+	return rgb*65536 + rgb*256 + rgb;
+}
+
+extern "C"
+{
+
+__global__ void translate(int max_idx, float3 *x, float3 v)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	x[idx] += v;
+}
+
+__global__ void rotate(int max_idx, float3 *x, float phi, float3 axis)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	x[idx] = rotate(x[idx], phi, axis);
+}
+
+__global__ void intersect_triangle_mesh(int max_idx, float3 *x, float3 *p, int n, float3 *mesh, int *pixel)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	bool hit = false;
+
+	float distance, min_distance;
+	float3 intersection, min_intersection;
+
+	int i;
+	for (i=0; i < n; i++)
+	{
+		if (intersect_triangle(x[idx], p[idx], mesh+3*i, intersection))
+		{
+			if (!hit)
+			{
+				pixel[idx] = get_color(p[idx], mesh+3*i);
+
+				min_distance = norm(intersection-x[idx]);
+				min_intersection = intersection;
+				hit = true;
+				continue;
+			}
+
+			distance = norm(intersection-x[idx]);
+
+			if (distance < min_distance)
+			{
+				pixel[idx] = get_color(p[idx], mesh+3*i);
+
+				min_distance = distance;
+				min_intersection = intersection;
+			}
+		}
+	}
+	
+	if (!hit)
+		pixel[idx] = 0;
+	
+}
+
+} // extern "c"
diff --git a/src/linalg.h b/src/linalg.h
new file mode 100644
index 0000000..593aa9d
--- /dev/null
+++ b/src/linalg.h
@@ -0,0 +1,116 @@
+#ifndef __LINALG_H__
+#define __LINALG_H__
+
+__device__ __host__ float3 operator- (const float3 &a)
+{
+	return make_float3(-a.x, -a.y, -a.z);
+}
+
+__device__ __host__ float3 operator+ (const float3 &a, const float3 &b)
+{
+	return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
+}
+
+__device__ __host__ void operator+= (float3 &a, const float3 &b)
+{
+	a.x += b.x;
+	a.y += b.y;
+	a.z += b.z;
+}
+
+__device__ __host__ float3 operator- (const float3 &a, const float3 &b)
+{
+	return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
+}
+
+__device__ __host__ void operator-= (float3 &a, const float3 &b)
+{
+	a.x -= b.x;
+	a.y -= b.y;
+	a.z -= b.z;
+}
+
+__device__ __host__ float3 operator+ (const float3 &a, const float &c)
+{
+	return make_float3(a.x+c, a.y+c, a.z+c);
+}
+
+__device__ __host__ void operator+= (float3 &a, const float &c)
+{
+	a.x += c;
+	a.y += c;
+	a.z += c;
+}
+
+__device__ __host__ float3 operator+ (const float &c, const float3 &a)
+{
+	return make_float3(c+a.x, c+a.y, c+a.z);
+}
+
+__device__ __host__ float3 operator- (const float3 &a, const float &c)
+{
+	return make_float3(a.x-c, a.y-c, a.z-c);
+}
+
+__device__ __host__ void operator-= (float3 &a, const float &c)
+{
+	a.x -= c;
+	a.y -= c;
+	a.z -= c;
+}
+
+__device__ __host__ float3 operator- (const float &c, const float3& a)
+{
+	return make_float3(c-a.x, c-a.y, c-a.z);
+}
+
+__device__ __host__ float3 operator* (const float3 &a, const float &c)
+{
+	return make_float3(a.x*c, a.y*c, a.z*c);
+}
+
+__device__ __host__ void operator*= (float3 &a, const float &c)
+{
+	a.x *= c;
+	a.y *= c;
+	a.z *= c;
+}
+
+__device__ __host__ float3 operator* (const float &c, const float3& a)
+{
+	return make_float3(c*a.x, c*a.y, c*a.z);
+}
+
+__device__ __host__ float3 operator/ (const float3 &a, const float &c)
+{
+	return make_float3(a.x/c, a.y/c, a.z/c);
+}
+
+__device__ __host__ void operator/= (float3 &a, const float &c)
+{
+	a.x /= c;
+	a.y /= c;
+	a.z /= c;
+}
+
+__device__ __host__ float3 operator/ (const float &c, const float3 &a)
+{
+	return make_float3(c/a.x, c/a.y, c/a.z);
+}
+
+__device__ __host__ float dot(const float3 &a, const float3 &b)
+{
+	return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+__device__ __host__ float3 cross(const float3 &a, const float3 &b)
+{
+	return make_float3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
+}
+
+__device__ __host__ float norm(const float3 &a)
+{
+	return sqrtf(dot(a,a));
+}
+
+#endif
diff --git a/src/matrix.h b/src/matrix.h
new file mode 100644
index 0000000..5c72344
--- /dev/null
+++ b/src/matrix.h
@@ -0,0 +1,225 @@
+#ifndef __MATRIX_H__
+#define __MATRIX_H__
+
+struct Matrix
+{
+	float a00, a01, a02, a10, a11, a12, a20, a21, a22;
+};
+
+__device__ __host__ Matrix make_matrix(float a00, float a01, float a02,
+				       float a10, float a11, float a12,
+				       float a20, float a21, float a22)
+{
+	Matrix m = {a00, a01, a02, a10, a11, a12, a20, a21, a22};
+	return m;
+}
+
+__device__ __host__ Matrix make_matrix(const float3 &u1, const float3 &u2, const float3 &u3)
+{
+	Matrix m = {u1.x, u2.x, u3.x, u1.y, u2.y, u3.y, u1.z, u2.z, u3.z};
+	return m;
+}
+
+__device__ __host__ Matrix operator- (const Matrix &m)
+{
+	return make_matrix(-m.a00, -m.a01, -m.a02,
+			   -m.a10, -m.a11, -m.a12,
+			   -m.a20, -m.a21, -m.a22);
+}
+
+__device__ __host__ float3 operator* (const Matrix &m, const float3 &a)
+{
+	return make_float3(m.a00*a.x + m.a01*a.y + m.a02*a.z,
+			   m.a10*a.x + m.a11*a.y + m.a12*a.z,
+			   m.a20*a.x + m.a21*a.y + m.a22*a.z);
+}
+
+__device__ __host__ Matrix operator+ (const Matrix &m, const Matrix &n)
+{
+	return make_matrix(m.a00+n.a00, m.a01+n.a01, m.a02+n.a02,
+			   m.a10+n.a10, m.a11+n.a11, m.a12+n.a12,
+			   m.a20+n.a20, m.a21+n.a21, m.a22+n.a22);
+}
+
+__device__ __host__ void operator+= (Matrix &m, const Matrix &n)
+{
+	m.a00 += n.a00;
+	m.a01 += n.a01;
+	m.a02 += n.a02;
+	m.a10 += n.a10;
+	m.a11 += n.a11;
+	m.a12 += n.a12;
+	m.a20 += n.a20;
+	m.a21 += n.a21;
+	m.a22 += n.a22;
+}
+
+__device__ __host__ Matrix operator- (const Matrix &m, const Matrix &n)
+{
+	return make_matrix(m.a00-n.a00, m.a01-n.a01, m.a02-n.a02,
+			   m.a10-n.a10, m.a11-n.a11, m.a12-n.a12,
+			   m.a20-n.a20, m.a21-n.a21, m.a22-n.a22);
+}
+
+__device__ __host__ void operator-= (Matrix &m, const Matrix& n)
+{
+	m.a00 -= n.a00;
+	m.a01 -= n.a01;
+	m.a02 -= n.a02;
+	m.a10 -= n.a10;
+	m.a11 -= n.a11;
+	m.a12 -= n.a12;
+	m.a20 -= n.a20;
+	m.a21 -= n.a21;
+	m.a22 -= n.a22;
+}
+
+__device__ __host__ Matrix operator* (const Matrix &m, const Matrix &n)
+{
+	return make_matrix(m.a00*n.a00 + m.a01*n.a10 + m.a02*n.a20,
+			   m.a00*n.a01 + m.a01*n.a11 + m.a02*n.a21,
+			   m.a00*n.a02 + m.a01*n.a12 + m.a02*n.a22,
+			   m.a10*n.a00 + m.a11*n.a10 + m.a12*n.a20,
+			   m.a10*n.a01 + m.a11*n.a11 + m.a12*n.a21,
+			   m.a10*n.a02 + m.a11*n.a12 + m.a12*n.a22,
+			   m.a20*n.a00 + m.a21*n.a10 + m.a22*n.a20,
+			   m.a20*n.a01 + m.a21*n.a11 + m.a22*n.a21,
+			   m.a20*n.a02 + m.a21*n.a12 + m.a22*n.a22);
+}
+
+__device__ __host__ Matrix operator+ (const Matrix &m, const float &c)
+{
+	return make_matrix(m.a00+c, m.a01+c, m.a02+c,
+			   m.a10+c, m.a11+c, m.a12+c,
+			   m.a20+c, m.a21+c, m.a22+c);
+}
+
+__device__ __host__ void operator+= (Matrix &m, const float &c)
+{
+	m.a00 += c;
+	m.a01 += c;
+	m.a02 += c;
+	m.a10 += c;
+	m.a11 += c;
+	m.a12 += c;
+	m.a20 += c;
+	m.a21 += c;
+	m.a22 += c;
+}
+
+__device__ __host__ Matrix operator+ (const float &c, const Matrix &m)
+{
+	return make_matrix(c+m.a00, c+m.a01, c+m.a02,
+			   c+m.a10, c+m.a11, c+m.a12,
+			   c+m.a20, c+m.a21, c+m.a22);
+}
+
+__device__ __host__ Matrix operator- (const Matrix &m, const float &c)
+{
+	return make_matrix(m.a00-c, m.a01-c, m.a02-c,
+			   m.a10-c, m.a11-c, m.a12-c,
+			   m.a20-c, m.a21-c, m.a22-c);
+}
+
+__device__ __host__ void operator-= (Matrix &m, const float &c)
+{
+	m.a00 -= c;
+	m.a01 -= c;
+	m.a02 -= c;
+	m.a10 -= c;
+	m.a11 -= c;
+	m.a12 -= c;
+	m.a20 -= c;
+	m.a21 -= c;
+	m.a22 -= c;
+}
+
+__device__ __host__ Matrix operator- (const float &c, const Matrix &m)
+{
+	return make_matrix(c-m.a00, c-m.a01, c-m.a02,
+			   c-m.a10, c-m.a11, c-m.a12,
+			   c-m.a20, c-m.a21, c-m.a22);
+}
+
+__device__ __host__ Matrix operator* (const Matrix &m, const float &c)
+{
+	return make_matrix(m.a00*c, m.a01*c, m.a02*c,
+			   m.a10*c, m.a11*c, m.a12*c,
+			   m.a20*c, m.a21*c, m.a22*c);
+}
+
+__device__ __host__ void operator*= (Matrix &m, const float &c)
+{
+	m.a00 *= c;
+	m.a01 *= c;
+	m.a02 *= c;
+	m.a10 *= c;
+	m.a11 *= c;
+	m.a12 *= c;
+	m.a20 *= c;
+	m.a21 *= c;
+	m.a22 *= c;
+}
+
+__device__ __host__ Matrix operator* (const float &c, const Matrix &m)
+{
+	return make_matrix(c*m.a00, c*m.a01, c*m.a02,
+			   c*m.a10, c*m.a11, c*m.a12,
+			   c*m.a20, c*m.a21, c*m.a22);
+}
+
+__device__ __host__ Matrix operator/ (const Matrix &m, const float &c)
+{
+	return make_matrix(m.a00/c, m.a01/c, m.a02/c,
+			   m.a10/c, m.a11/c, m.a12/c,
+			   m.a20/c, m.a21/c, m.a22/c);
+}
+
+__device__ __host__ void operator/= (Matrix &m, const float &c)
+{
+	m.a00 /= c;
+	m.a01 /= c;
+	m.a02 /= c;
+	m.a10 /= c;
+	m.a11 /= c;
+	m.a12 /= c;
+	m.a20 /= c;
+	m.a21 /= c;
+	m.a22 /= c;
+}
+
+__device__ __host__ Matrix operator/ (const float &c, const Matrix &m)
+{
+	return make_matrix(c/m.a00, c/m.a01, c/m.a02,
+			   c/m.a10, c/m.a11, c/m.a12,
+			   c/m.a20, c/m.a21, c/m.a22);
+}
+
+__device__ __host__ float det(const Matrix &m)
+{
+	return m.a00*(m.a11*m.a22 - m.a12*m.a21) - \
+		m.a10*(m.a01*m.a22 - m.a02*m.a21) + \
+		m.a20*(m.a01*m.a12 - m.a02*m.a11);
+}
+
+__device__ __host__ Matrix inv(const Matrix &m)
+{
+	return make_matrix(m.a11*m.a22 - m.a12*m.a21,
+			   m.a02*m.a21 - m.a01*m.a22,
+			   m.a01*m.a12 - m.a02*m.a11,
+			   m.a12*m.a20 - m.a10*m.a22,
+			   m.a00*m.a22 - m.a02*m.a20,
+			   m.a02*m.a10 - m.a00*m.a12,
+			   m.a10*m.a21 - m.a11*m.a20,
+			   m.a01*m.a20 - m.a00*m.a21,
+			   m.a00*m.a11 - m.a01*m.a10)/det(m);
+}
+
+__device__ __host__ Matrix outer(const float3 &a, const float3 &b)
+{
+	return make_matrix(a.x*b.x, a.x*b.y, a.x*b.z,
+			   a.y*b.x, a.y*b.y, a.y*b.z,
+			   a.z*b.x, a.z*b.y, a.z*b.z);
+}
+
+#endif
diff --git a/src/rotate.h b/src/rotate.h
new file mode 100644
index 0000000..52d6d6a
--- /dev/null
+++ b/src/rotate.h
@@ -0,0 +1,24 @@
+#ifndef __ROTATE_H__
+#define __ROTATE_H__
+
+__device__ const Matrix IDENTITY_MATRIX = {1,0,0,0,1,0,0,0,1};
+__device__ const Matrix ZERO_MATRIX = {0,0,0,0,0,0,0,0,0};
+
+__device__ __host__ Matrix make_rotation_matrix(float phi, const float3 &n)
+{
+	/* rotate points counterclockwise, when looking towards +infinity,
+	   through an angle `phi` about the axis `n`. */
+
+	float cos_phi = cosf(phi);
+	float sin_phi = sinf(phi);
+
+	return IDENTITY_MATRIX*cos_phi + (1-cos_phi)*outer(n,n) +
+		sin_phi*make_matrix(0,n.z,-n.y,-n.z,0,n.x,n.y,-n.x,0);
+}
+
+__device__ __host__ float3 rotate(const float3 &a, float phi, const float3 &n)
+{
+	return make_rotation_matrix(phi,n)*a;
+}
+
+#endif
diff --git a/stl.py b/stl.py
index b300554..076c739 100644
--- a/stl.py
+++ b/stl.py
@@ -1,15 +1,30 @@
 import numpy as np
+import struct
 
 def pull_vertices(filename):
     f = open(filename)
 
-    vertices = []
+    vertex = []
     for line in f:
         if not line.strip().startswith('vertex'):
             continue
-        vertices.append([float(s) for s in line.strip().split()[1:]])
+        vertex.append([float(s) for s in line.strip().split()[1:]])
 
-    return np.array(vertices)
+    f.close()
+    return np.array(vertex)
 
-if __name__ == '__main__':
-    print pull_vertices('models/MiniFig.STL')
+def pull_vertices_binary(filename):
+    f = open(filename)
+
+    f.read(80)
+    triangles = struct.unpack('<I', f.read(4))[0]
+
+    vertex = []
+    for i in range(triangles):
+        f.read(12)
+        for j in range(3):
+            vertex.append(struct.unpack('<fff', f.read(12)))
+        f.read(2)
+
+    f.close()
+    return np.array(vertex)
diff --git a/test.py b/test.py
new file mode 100644
index 0000000..014c2db
--- /dev/null
+++ b/test.py
@@ -0,0 +1,88 @@
+import time
+from stl import *
+import numpy as np
+from pycuda import autoinit
+from pycuda.compiler import SourceModule
+import pycuda.driver as cuda
+from pycuda import gpuarray
+
+def array2float3(arr):
+    if len(arr.shape) != 2 or arr.shape[-1] != 3:
+        raise Exception('shape mismatch')
+
+    x = np.empty(arr.shape[0], dtype=gpuarray.vec.float3)
+    x['x'] = arr[:,0]
+    x['y'] = arr[:,1]
+    x['z'] = arr[:,2]
+
+    return x
+
+print 'device %s' % autoinit.device.name()
+
+source = open('src/linalg.h').read() + open('src/matrix.h').read() + \
+    open('src/rotate.h').read() + open('src/intersect.cu').read()
+
+mod = SourceModule(source, no_extern_c=True, arch='sm_13')
+intersect = mod.get_function('intersect_triangle_mesh')
+rotate = mod.get_function('rotate')
+translate = mod.get_function('translate')
+
+mesh = array2float3(pull_vertices_binary('models/tie_interceptor6.stl'))
+
+import pygame
+size = width, height = 800, 600
+screen = pygame.display.set_mode(size)
+
+film_size = (0.035, 0.024)
+focal_length = 0.05
+
+grid = []
+for x in np.linspace(-film_size[0]/2, film_size[0]/2, width):
+    for z in np.linspace(-film_size[1]/2, film_size[1]/2, height):
+        grid.append((x,0,z))
+grid = np.array(grid)
+grid += (0,focal_length,0)
+grid += (0,300,0)
+
+x = array2float3(grid)
+p = array2float3(((0,300,0)-grid))
+
+x_gpu = cuda.mem_alloc(x.nbytes)
+cuda.memcpy_htod(x_gpu,x)
+
+p_gpu = cuda.mem_alloc(p.nbytes)
+cuda.memcpy_htod(p_gpu,p)
+
+mesh_gpu = cuda.mem_alloc(mesh.nbytes)
+cuda.memcpy_htod(mesh_gpu,mesh)
+
+pixel = np.empty(size, dtype=np.int32).flatten()
+pixel_gpu = cuda.mem_alloc(pixel.nbytes)
+cuda.memcpy_htod(pixel_gpu,pixel)
+
+rotate(np.int32(mesh.size), mesh_gpu, np.float32(-np.pi/2), gpuarray.vec.make_float3(1,0,0), block=(256,1,1), grid=(mesh.size//256+1,1))
+
+translate(np.int32(mesh.size), mesh_gpu, gpuarray.vec.make_float3(0,30,0), block=(256,1,1), grid=(mesh.size//256+1,1))
+
+t0 = time.time()
+for i in range(100):
+    rotate(np.int32(x.size), x_gpu, np.float32(np.pi/50), gpuarray.vec.make_float3(0,0,1), block=(256,1,1), grid=(width*height//256+1,1))
+
+    rotate(np.int32(p.size), p_gpu, np.float32(np.pi/50), gpuarray.vec.make_float3(0,0,1), block=(256,1,1), grid=(width*height//256+1,1))
+
+    intersect(np.int32(x.size), x_gpu, p_gpu, np.int32(mesh.size//3), mesh_gpu, pixel_gpu, block=(256,1,1), grid=(width*height//256+1,1))
+
+    cuda.Context.synchronize()
+
+    cuda.memcpy_dtoh(pixel, pixel_gpu)
+    pygame.surfarray.blit_array(screen, pixel.reshape(size))
+    pygame.display.flip()
+
+
+elapsed = time.time() - t0
+
+print '%i triangles, %i photons, %f sec; (%f photons/s)' % \
+    (mesh.size//3, pixel.size, elapsed, pixel.size/elapsed)
+
+
+raw_input('press enter to exit')
-- 
cgit