tie fighter

4 files changed, 483 insertions, 0 deletions

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