added some more documentation and a more accurate miniature version of lbne

2 files changed, 193 insertions, 159 deletions

diff --git a/src/intersect.cu b/src/intersect.h
index 1402aa1..b984612 100644
--- a/src/intersect.cu
+++ b/src/intersect.h
@@ -5,22 +5,6 @@
 #include "matrix.h"
 #include "rotate.h"
 
-/* flattened triangle mesh */
-texture<float4, 1, cudaReadModeElementType> mesh;
-
-/* lower/upper bounds for the bounding box associated with each node/leaf */
-texture<float4, 1, cudaReadModeElementType> upper_bound_arr;
-texture<float4, 1, cudaReadModeElementType> lower_bound_arr;
-
-/* map to child nodes/triangles and the number of child nodes/triangles */
-texture<uint, 1, cudaReadModeElementType> child_map_arr;
-texture<uint, 1, cudaReadModeElementType> child_len_arr;
-
-__device__ float3 make_float3(const float4 &a)
-{
-	return make_float3(a.x, a.y, a.z);
-}
-
 /* Test the intersection between a ray starting from `origin` traveling in the
    direction `direction` and a triangle defined by the vertices `v0`, `v1`, and
    `v2`. If the ray intersects the triangle, set `distance` to the distance
@@ -155,146 +139,3 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con
 
 	return true;
 }
-
-/* Test the intersection between a ray starting at `origin` traveling in the
-   direction `direction` and the bounding box around node `i`. If the ray
-   intersects the bounding box return true, else return false. */
-__device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i)
-{
-	float3 lower_bound = make_float3(tex1Dfetch(lower_bound_arr, i));
-	float3 upper_bound = make_float3(tex1Dfetch(upper_bound_arr, i));
-
-	return intersect_box(origin, direction, lower_bound, upper_bound);
-}
-
-extern "C"
-{
-
-__global__ void translate(int max_idx, float3 *pt, float3 v)
-{
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
-
-	if (idx > max_idx)
-		return;
-
-	pt[idx] += v;
-}
-
-__global__ void rotate(int max_idx, float3 *pt, float phi, float3 axis)
-{
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
-
-	if (idx > max_idx)
-		return;
-
-	pt[idx] = rotate(pt[idx], phi, axis);
-}
-
-__global__ void intersect_mesh(int max_idx, float3 *origin_arr, float3 *direction_arr, int first_leaf, int *state_arr, int *pixel_arr)
-{
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
-
-	if (idx > max_idx)
-		return;
-
-	float3 origin = origin_arr[idx];
-	float3 direction = direction_arr[idx];
-	direction /= norm(direction);
-
-	int *pixel = pixel_arr+idx;
-	int *state = state_arr+idx;
-
-	bool hit = false;
-
-	float distance;
-	float min_distance;
-	
-	if (!intersect_node(origin, direction, 0))
-	{
-		*pixel = 0;
-		return;
-	}
-
-	int stack[100];
-
-	int *head = &stack[0];
-	int *node = &stack[1];
-	*node = 0;
-
-	int i;
-
-	bool show_leafs = false;
-
-	do
-	{
-		int child_map = tex1Dfetch(child_map_arr, *node);
-		int child_len = tex1Dfetch(child_len_arr, *node);
-
-		if (*node < first_leaf)
-		{
-			for (i=0; i < child_len; i++)
-				if (intersect_node(origin, direction, child_map+i))
-					*node++ = child_map+i;
-
-			if (node == head)
-				break;
-			
-			node--;
-
-		}
-		else if (show_leafs)
-		{
-			hit = true;
-			*pixel = 255;
-			node--;
-		}
-		else // node is a leaf
-		{
-			for (i=0; i < child_len; i++)
-			{
-				int mesh_idx = 3*(child_map + i);
-
-				float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx));
-				float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1));
-				float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2));
-				
-				if (intersect_triangle(origin, direction, v0, v1, v2, distance))
-				{
-					if (!hit)
-					{
-						*pixel = get_color(direction, v0, v1, v2);
-						*state = child_map + i;
-
-						min_distance = distance;
-
-						hit = true;
-						continue;
-					}
-					
-					if (distance < min_distance)
-					{
-						*pixel = get_color(direction, v0, v1, v2);
-						*state = child_map + i;
-						
-						min_distance = distance;
-					}
-				}
-
-			} // triangle loop
-
-			node--;
-
-		} // node is a leaf
-
-	} // while loop
-	while (node != head);
-
-	if (!hit)
-	{
-		*state = -1;
-		*pixel = 0;
-	}
-
-} // intersect mesh
-
-} // extern "c"
diff --git a/src/kernel.cu b/src/kernel.cu
new file mode 100644
index 0000000..ed53032
--- /dev/null
+++ b/src/kernel.cu
@@ -0,0 +1,193 @@
+//-*-c-*-
+#include <math_constants.h>
+
+#include "linalg.h"
+#include "matrix.h"
+#include "rotate.h"
+#include "intersect.h"
+
+#define STACK_SIZE 100
+
+/* flattened triangle mesh */
+texture<float4, 1, cudaReadModeElementType> mesh;
+
+/* lower/upper bounds for the bounding box associated with each node/leaf */
+texture<float4, 1, cudaReadModeElementType> upper_bounds;
+texture<float4, 1, cudaReadModeElementType> lower_bounds;
+
+/* map to child nodes/triangles and the number of child nodes/triangles */
+texture<uint, 1, cudaReadModeElementType> child_map_arr;
+texture<uint, 1, cudaReadModeElementType> child_len_arr;
+
+__device__ float3 make_float3(const float4 &a)
+{
+	return make_float3(a.x, a.y, a.z);
+}
+
+/* Test the intersection between a ray starting at `origin` traveling in the
+   direction `direction` and the bounding box around node `i`. If the ray
+   intersects the bounding box return true, else return false. */
+__device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i)
+{
+	float3 lower_bound = make_float3(tex1Dfetch(lower_bounds, i));
+	float3 upper_bound = make_float3(tex1Dfetch(upper_bounds, i));
+
+	return intersect_box(origin, direction, lower_bound, upper_bound);
+}
+
+/* Find the intersection between a ray starting at `origin` traveling in the
+   direction `direction` and the global mesh texture. If the ray intersects
+   the texture return the index of the triangle which the ray intersected,
+   else return -1. */
+__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int first_leaf)
+{
+	int triangle_idx = -1;
+
+	float distance;
+	float min_distance;
+
+	if (!intersect_node(origin, direction, 0))
+		return -1;
+
+	int stack[STACK_SIZE];
+
+	int *head = &stack[0];
+	int *node = &stack[1];
+	int *tail = &stack[STACK_SIZE-1];
+	*node = 0;
+
+	int i;
+
+	do
+	{
+		int child_map = tex1Dfetch(child_map_arr, *node);
+		int child_len = tex1Dfetch(child_len_arr, *node);
+
+		if (*node < first_leaf)
+		{
+			for (i=0; i < child_len; i++)
+				if (intersect_node(origin, direction, child_map+i))
+					*node++ = child_map+i;
+
+			if (node == head)
+				break;
+
+			node--;
+		}
+		else // node is a leaf
+		{
+			for (i=0; i < child_len; i++)
+			{
+				int mesh_idx = 3*(child_map + i);
+
+				float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx));
+				float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1));
+				float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2));
+
+				if (intersect_triangle(origin, direction, v0, v1, v2, distance))
+				{
+					if (triangle_idx == -1)
+					{
+						triangle_idx = child_map + i;
+						min_distance = distance;
+						continue;
+					}
+
+					if (distance < min_distance)
+					{
+						triangle_idx = child_map + i;
+						min_distance = distance;
+					}
+				}
+			} // triangle loop
+
+			node--;
+
+		} // node is a leaf
+
+	} // while loop
+	while (node != head);
+
+	return triangle_idx;
+}
+
+extern "C"
+{
+
+/* Translate `points` by the vector `v` */
+__global__ void translate(int max_idx, float3 *points, float3 v)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	*(points+idx) += v;
+}
+
+/* Rotate `points` through an angle `phi` counter-clockwise about the
+   axis `axis` (when looking towards +infinity). */
+__global__ void rotate(int max_idx, float3 *points, float phi, float3 axis)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	*(points+idx) = rotate(*(points+idx), phi, axis);
+}
+
+/* Trace the rays starting at `origins` traveling in the direction `directions`
+   to their intersection with the global mesh. If the ray intersects the mesh
+   set the pixel associated with the ray to a 32 bit color whose brightness is
+   determined by the cosine of the angle between the ray and the normal of the
+   triangle it intersected, else set the pixel to 0. */
+__global__ void ray_trace(int max_idx, float3 *origins, float3 *directions, int first_leaf, int *pixels)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	float3 origin = *(origins+idx);
+	float3 direction = *(directions+idx);
+	direction /= norm(direction);
+
+	int intersection_idx = intersect_mesh(origin, direction, first_leaf);
+
+	if (intersection_idx == -1)
+	{
+		*(pixels+idx) = 0;
+	}
+	else
+	{
+		int mesh_idx = 3*intersection_idx;
+
+		float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx));
+		float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1));
+		float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2));
+
+		*(pixels+idx) = get_color(direction, v0, v1, v2);
+	}
+} // ray_trace
+
+/* Propagate the photons starting at `origins` traveling in the direction
+   `directions` to their intersection with the global mesh. If the ray
+   intersects the mesh set the hit_solid array value associated with the
+   photon to the triangle index of the triangle the photon intersected, else
+   set the hit_solid array value to -1. */
+__global__ void propagate(int max_idx, float3 *origins, float3 *directions, int first_leaf, int *hit_solids)
+{
+	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (idx > max_idx)
+		return;
+
+	float3 origin = *(origins+idx);
+	float3 direction = *(directions+idx);
+	direction /= norm(direction);
+
+	*(hit_solids+idx) = intersect_mesh(origin, direction, first_leaf);
+} // propagate
+
+} // extern "c"