1 files changed, 72 insertions, 0 deletions
diff --git a/intersect.cu b/intersect.cu
new file mode 100644
index 0000000..d7fdaa5
--- /dev/null
+++ b/intersect.cu
@@ -0,0 +1,72 @@
+//-*-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"