pulled a lot of the photon propagation code out of src/kernel.cu into src/photon.h so that photon propagation by propagate() in kernel.cu and the hybrid monte carlo ray tracing use the same code. instead of a single state, photons now carry the history of the processes they've undergone. this history is stored as a bitmask; see src/photon.h. start_node and first_node of the mesh are now stored as global variables in mesh.h instead of being passed to kernel functions.

4 files changed, 548 insertions, 580 deletions

diff --git a/src/daq.cu b/src/daq.cu
index a7f5120..07e2d5e 100644
--- a/src/daq.cu
+++ b/src/daq.cu
@@ -16,7 +16,8 @@ __device__ float sortable_int_to_float(unsigned int i)
   //return __int_as_float(i ^ mask);
 }
 
-extern "C" {
+extern "C"
+{
 
 __global__ void reset_earliest_time_int(float maxtime, int ntime_ints, unsigned int *time_ints)
 {
@@ -27,8 +28,10 @@ __global__ void reset_earliest_time_int(float maxtime, int ntime_ints, unsigned
 	}
 }
 
-__global__ void run_daq(curandState *s, int detection_state, float time_rms,
-			int nphotons, float *photon_times, int *photon_states,
+__global__ void run_daq(curandState *s, unsigned int detection_state,
+			float time_rms,
+			int nphotons, float *photon_times,
+			unsigned int *photon_histories,
 			int *last_hit_triangles, int *solid_map,
 			int nsolids, unsigned int *earliest_time_int)
 {
@@ -39,37 +42,26 @@ __global__ void run_daq(curandState *s, int detection_state, float time_rms,
 	{
 		curandState rng = s[id];
 
-
-
 		int triangle_id = last_hit_triangles[id];
 		
 		if (triangle_id > -1)
 		{
 			int solid_id = solid_map[triangle_id];
-			int state = photon_states[id];
+			int history = photon_histories[id];
 			float time = photon_times[id] + curand_normal(&rng) * time_rms;
 			unsigned int time_int = float_to_sortable_int(time);
 
-
-
-			if (solid_id < nsolids && state == detection_state) {
+			if (solid_id < nsolids && (history & detection_state)) 
+			{
 			    atomicMin(earliest_time_int + solid_id, time_int);
 			}
 
-
 		}
 
-
-
 		s[id] = rng;
 
-
-
 	}
 
-
-
-
 }
 
 __global__ void convert_sortable_int_to_float(int n, 
diff --git a/src/kernel.cu b/src/kernel.cu
index 44a6c03..73800a3 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -5,364 +5,71 @@
 #include "linalg.h"
 #include "matrix.h"
 #include "rotate.h"
-#include "intersect.h"
-#include "materials.h"
-#include "physical_constants.h"
-#include "random.h"
-
-#define STACK_SIZE 500
+#include "mesh.h"
+#include "photon.h"
 
 #define RED_WAVELENGTH 685
 #define BLUE_WAVELENGTH 465
 #define GREEN_WAVELENGTH 545
 
-enum
-{
-	DEBUG = -2,
-	INIT = -1,
-	NO_HIT,
-	BULK_ABSORB,
-	SURFACE_DETECT,
-	SURFACE_ABSORB
-};
-
-/* flattened triangle mesh */
-__device__ float3 *vertices;
-__device__ uint4 *triangles;
-
-/* 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<unsigned int, 1, cudaReadModeElementType> node_map;
-texture<unsigned int, 1, cudaReadModeElementType> node_length;
-
-__device__ float3 make_float3(const float4 &a)
-{
-	return make_float3(a.x, a.y, a.z);
-}
-
-__device__ int convert(int c)
-{
-	if (c & 0x80)
-		return (0xFFFFFF00 | c);
-	else
-		return c;
-}
-
-/* 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 start_node, const int first_node, float &min_distance, int last_hit_triangle = -1)
-{
-	int triangle_index = -1;
-
-	float distance;
-
-	if (!intersect_node(origin, direction, start_node))
-		return -1;
-
-	int stack[STACK_SIZE];
-
-	int *head = &stack[0];
-	int *node = &stack[1];
-	int *tail = &stack[STACK_SIZE-1];
-	*node = start_node;
-
-	int i;
-
-	do
-	{
-		int first_child = tex1Dfetch(node_map, *node);
-		int child_count = tex1Dfetch(node_length, *node);
-
-		while (*node >= first_node && child_count == 1)
-		{
-			*node = first_child;
-			first_child = tex1Dfetch(node_map, *node);
-			child_count = tex1Dfetch(node_length, *node);
-		}
-		
-		if (*node >= first_node)
-		{
-			for (i=0; i < child_count; i++)
-			{
-				if (intersect_node(origin, direction, first_child+i))
-				{
-					*node = first_child+i;
-					node++;
-				}
-			}
-
-			node--;
-		}
-		else // node is a leaf
-		{
-			for (i=0; i < child_count; i++)
-			{
-				if (last_hit_triangle == first_child+i)
-					continue;
-
-				uint4 triangle_data = triangles[first_child+i];
-
-				float3 v0 = vertices[triangle_data.x];
-				float3 v1 = vertices[triangle_data.y];
-				float3 v2 = vertices[triangle_data.z];
-
-				if (intersect_triangle(origin, direction, v0, v1, v2, distance))
-				{
-					if (triangle_index == -1)
-					{
-						triangle_index = first_child + i;
-						min_distance = distance;
-						continue;
-					}
-
-					if (distance < min_distance)
-					{
-						triangle_index = first_child + i;
-						min_distance = distance;
-					}
-				}
-			} // triangle loop
-
-			node--;
-
-		} // node is a leaf
-
-	} // while loop
-	while (node != head);
-
-	return triangle_index;
-}
-
 __device__ void myAtomicAdd(float *addr, float data)
 {
 	while (data)
 		data = atomicExch(addr, data+atomicExch(addr, 0.0f));
 }
 
-__device__ int to_diffuse(curandState &rng, float3 position, float3 direction, float wavelength, float3 polarization, int start_node, int first_node, int max_steps)
+__device__ int to_diffuse(Photon p, int max_steps)
 {
-	int last_hit_triangle = -1;
+	// note that p is NOT passed by reference; this is intentional
+
+	p.last_hit_triangle = -1;
+
+	State s;
 
 	int steps = 0;
 	while (steps < max_steps)
 	{
 		steps++;
 
-		float distance;
-
-	        last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
-
-		if (last_hit_triangle == -1)
-			return last_hit_triangle;
-
-		uint4 triangle_data = triangles[last_hit_triangle];
-
-		float3 v0 = vertices[triangle_data.x];
-		float3 v1 = vertices[triangle_data.y];
-		float3 v2 = vertices[triangle_data.z];
+		int command;
 
-		int material_in_index = convert(0xFF & (triangle_data.w >> 24));
-		int material_out_index = convert(0xFF & (triangle_data.w >> 16));
-		int surface_index = convert(0xFF & (triangle_data.w >> 8));
+		command = fill_state(s, p);
 
-		float3 surface_normal = cross(v1-v0,v2-v1);
-		surface_normal /= norm(surface_normal);
-		
-		Material material1, material2;
-		if (dot(surface_normal,-direction) > 0.0f)
-		{
-			// outside to inside
-			material1 = materials[material_out_index];
-			material2 = materials[material_in_index];
-		}
-		else
-		{
-			// inside to outside
-			material1 = materials[material_in_index];
-			material2 = materials[material_out_index];
-			surface_normal = -surface_normal;
-		}
+		if (command == BREAK)
+			break;
 
-		float refractive_index1 = interp_property(wavelength, material1.refractive_index);
-		float refractive_index2 = interp_property(wavelength, material2.refractive_index);
+		command = propagate_to_boundary(p, s);
 
-		float absorption_length = interp_property(wavelength, material1.absorption_length);
-		float scattering_length = interp_property(wavelength, material1.scattering_length);
+		if (command == BREAK)
+			break;
 
-		float absorption_distance = -absorption_length*logf(curand_uniform(&rng));
-		float scattering_distance = -scattering_length*logf(curand_uniform(&rng));
+		if (command == CONTINUE)
+			continue;
 
-		if (absorption_distance <= scattering_distance)
-		{
-			if (absorption_distance <= distance)
-			{
-				return -1;
-			} // photon is absorbed in material1
-		}
-		else
+		if (s.surface_index != -1)
 		{
-			if (scattering_distance <= distance)
-			{
-				//time += scattering_distance/(SPEED_OF_LIGHT/refractive_index1);
-				position += scattering_distance*direction;
-
-				float theta,y;
-				while (true)
-				{
-					y = curand_uniform(&rng);
-					theta = uniform(&rng, 0, 2*PI);
-
-					if (y < powf(cosf(theta),2))
-						break;
-				}
-				
-				float phi = uniform(&rng, 0, 2*PI);
-
-				float3 b = cross(polarization, direction);
-				float3 c = polarization;
-
-				direction = rotate(direction, theta, b);
-				direction = rotate(direction, phi, c);
-				polarization = rotate(polarization, theta, b);
-				polarization = rotate(polarization, phi, c);
-
-				last_hit_triangle = -1;
-
-				continue;
-			} // photon is scattered in material1
-
-		} // if scattering_distance < absorption_distance
-
-		position += distance*direction;
-		//time += distance/(SPEED_OF_LIGHT/refractive_index1);
+			command = propagate_at_surface(p, s);
 
-		// p is normal to the plane of incidence
-		float3 p = cross(direction, surface_normal);
-		p /= norm(p);
+			if (p.history & REFLECT_DIFFUSE)
+				return p.last_hit_triangle;
 
-		float normal_coefficient = dot(polarization, p);
-		float normal_probability = normal_coefficient*normal_coefficient;
-
-		float incident_angle = acosf(dot(surface_normal,-direction));
-
-		if (surface_index != -1)
-		{
-			Surface surface = surfaces[surface_index];
-
-			float detect = interp_property(wavelength, surface.detect);
-			float absorb = interp_property(wavelength, surface.absorb);
-			float reflect_diffuse = interp_property(wavelength, surface.reflect_diffuse);
-			float reflect_specular = interp_property(wavelength, surface.reflect_specular);
-
-			// since the surface properties are interpolated linearly,
-			// we are guaranteed that they still sum to one.
-
-			float uniform_sample = curand_uniform(&rng);
-
-			if (uniform_sample < absorb)
-			{
-				// absorb
-				//state = SURFACE_ABSORB;
-				//break;
-				return -1;
-			}
-			else if (uniform_sample < absorb + detect)
-			{
-				// detect
-				//state = SURFACE_DETECT;
-				//break;
-				return -1;
-			}
-			else if (uniform_sample < absorb + detect + reflect_diffuse)
-			{
-				//state = DIFFUSE_HIT;
-				//break;
-				return last_hit_triangle;
-			}
-			else
-			{
-				// specularly reflect
-				direction = rotate(surface_normal, incident_angle, p);
-
-				// polarization ?
+			if (command == BREAK)
+				break;
 
+			if (command == CONTINUE)
 				continue;
-			}
-
-		} // surface
-
-		float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2);
-
-		if (curand_uniform(&rng) < normal_probability)
-		{
-
-			float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
-			float reflection_probability = reflection_coefficient*reflection_coefficient;
-
-			if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
-			{
-				direction = rotate(surface_normal, incident_angle, p);
-			}
-			else
-			{
-				direction = rotate(surface_normal, PI-refracted_angle, p);
-			}
-
-			polarization = p;
-
-			continue;
-		}  // photon polarization normal to plane of incidence
-		else
-		{
-			float reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
-			float reflection_probability = reflection_coefficient*reflection_coefficient;
-
-			if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
-			{
-				direction = rotate(surface_normal, incident_angle, p);
-			}
-			else
-			{
-				direction = rotate(surface_normal, PI-refracted_angle, p);
-			}
-
-			polarization = cross(p, direction);
-			polarization /= norm(polarization);
+		}
 
-			continue;
-		} // photon polarization parallel to surface
+		propagate_at_boundary(p, s);
 
-	} // while(nsteps < max_steps)
+	} // while (steps < max_steps)
 
 	return -1;
-
 } // to_diffuse
 
 extern "C"
 {
 
-__global__ void set_pointer(uint4 *triangle_ptr, float3 *vertex_ptr)
-{
-	triangles = triangle_ptr;
-	vertices = vertex_ptr;
-}
-
 /* Translate `points` by the vector `v` */
 __global__ void translate(int nthreads, float3 *points, float3 v)
 {
@@ -386,52 +93,61 @@ __global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
 	points[id] = rotate(points[id], phi, axis);
 }
 
-__global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, int start_node, int first_node, float3 *rgb_lookup, int runs, int max_steps)
+__global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float3 *rgb_lookup, int runs, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	curandState rng = rng_states[id];
-	float3 position = positions[id];
-	float3 direction = directions[id];
-	direction /= norm(direction);
-	float3 polarization = uniform_sphere(&rng);
+	Photon p;
+	p.rng = rng_states[id];
+	p.position = positions[id];
+	p.direction = directions[id];
+	p.direction /= norm(p.direction);
+	p.polarization = uniform_sphere(&p.rng);
+	p.time = 0.0f;
+	p.history = 0x0;
 
 	float distance;
 
-	int hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance);
+	int hit_triangle = intersect_mesh(p.position, p.direction, distance);
 
 	if (hit_triangle != id)
 		return;
 
 	// note triangles from built geometry mesh are in order
 
-	uint4 triangle_data = triangles[hit_triangle];
+	uint4 triangle_data = g_triangles[hit_triangle];
 
-	float3 v0 = vertices[triangle_data.x];
-	float3 v1 = vertices[triangle_data.y];
-	float3 v2 = vertices[triangle_data.z];
+	float3 v0 = g_vertices[triangle_data.x];
+	float3 v1 = g_vertices[triangle_data.y];
+	float3 v2 = g_vertices[triangle_data.z];
 
-	float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -direction);
+	float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -p.direction);
 
 	if (cos_theta < 0.0f)
-		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -direction);
+		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -p.direction);
 
 	for (int i=0; i < runs; i++)
 	{
-		hit_triangle = to_diffuse(rng, position, direction, RED_WAVELENGTH, polarization, start_node, first_node, max_steps);
+		p.wavelength = RED_WAVELENGTH;
+
+		hit_triangle = to_diffuse(p, max_steps);
 
 		if (hit_triangle != -1)
 			myAtomicAdd(&rgb_lookup[hit_triangle].x, cos_theta);
 
-		hit_triangle = to_diffuse(rng, position, direction, BLUE_WAVELENGTH, polarization, start_node, first_node, max_steps);
+		p.wavelength = BLUE_WAVELENGTH;
+
+		hit_triangle = to_diffuse(p, max_steps);
 
 		if (hit_triangle != -1)
 			myAtomicAdd(&rgb_lookup[hit_triangle].y, cos_theta);
 
-		hit_triangle = to_diffuse(rng, position, direction, GREEN_WAVELENGTH, polarization, start_node, first_node, max_steps);
+		p.wavelength = GREEN_WAVELENGTH;
+
+		hit_triangle = to_diffuse(p, max_steps);
 
 		if (hit_triangle != -1)
 			myAtomicAdd(&rgb_lookup[hit_triangle].z, cos_theta);
@@ -439,18 +155,21 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 
 } // build_rgb_lookup
 
-__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, int start_node, int first_node, float3 *rgb_lookup, int runs, int *pixels, int max_steps)
+__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float3 *rgb_lookup, int runs, int *pixels, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	curandState rng = rng_states[id];
-	float3 position = positions[id];
-	float3 direction = directions[id];
-	direction /= norm(direction);
-	float3 polarization = uniform_sphere(&rng);
+	Photon p;
+	p.rng = rng_states[id];
+	p.position = positions[id];
+	p.direction = directions[id];
+	p.direction /= norm(p.direction);
+	p.polarization = uniform_sphere(&p.rng);
+	p.time = 0.0f;
+	p.history = 0x0;
 
 	float3 rgb = make_float3(0.0, 0.0, 0.0);
 
@@ -458,17 +177,23 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 
 	for (int i=0; i < runs; i++)
 	{
-		hit_triangle = to_diffuse(rng, position, direction, RED_WAVELENGTH, polarization, start_node, first_node, max_steps);
+		p.wavelength = RED_WAVELENGTH;
+
+		hit_triangle = to_diffuse(p, max_steps);
 
 		if (hit_triangle != -1)
 			rgb.x += rgb_lookup[hit_triangle].x;
 
-		hit_triangle = to_diffuse(rng, position, direction, BLUE_WAVELENGTH, polarization, start_node, first_node, max_steps);
+		p.wavelength = BLUE_WAVELENGTH;
+
+		hit_triangle = to_diffuse(p, max_steps);
 
 		if (hit_triangle != -1)
 			rgb.y += rgb_lookup[hit_triangle].y;
 
-		hit_triangle = to_diffuse(rng, position, direction, GREEN_WAVELENGTH, polarization, start_node, first_node, max_steps);
+		p.wavelength = GREEN_WAVELENGTH;
+
+		hit_triangle = to_diffuse(p, max_steps);
 
 		if (hit_triangle != -1)
 			rgb.z += rgb_lookup[hit_triangle].z;
@@ -484,12 +209,13 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 
 } // render
 
-/* Trace the rays starting at `positions` 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 nthreads, float3 *positions, float3 *directions, int start_node, int first_node, int *pixels)
+/* Trace the rays starting at `positions` 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 nthreads, float3 *positions, float3 *directions, int *pixels)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
@@ -502,7 +228,7 @@ __global__ void ray_trace(int nthreads, float3 *positions, float3 *directions, i
 
 	float distance;
 
-	int triangle_index = intersect_mesh(position, direction, start_node, first_node, distance);
+	int triangle_index = intersect_mesh(position, direction, distance);
 
 	if (triangle_index == -1)
 	{
@@ -510,253 +236,84 @@ __global__ void ray_trace(int nthreads, float3 *positions, float3 *directions, i
 	}
 	else
 	{
-		uint4 triangle_data = triangles[triangle_index];
+		uint4 triangle_data = g_triangles[triangle_index];
 
-		float3 v0 = vertices[triangle_data.x];
-		float3 v1 = vertices[triangle_data.y];
-		float3 v2 = vertices[triangle_data.z];
+		float3 v0 = g_vertices[triangle_data.x];
+		float3 v1 = g_vertices[triangle_data.y];
+		float3 v2 = g_vertices[triangle_data.z];
 
 		pixels[id] = get_color(direction, v0, v1, v2, triangle_data.w);
 	}
 
 } // ray_trace
 
-__global__ void propagate(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, int *states, int *last_hit_triangles, int start_node, int first_node, int max_steps)
+__global__ void propagate(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, unsigned int *histories, int *last_hit_triangles, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	int state = states[id];
-
-	if (state != INIT)
+	Photon p;
+	p.rng = rng_states[id];
+	p.position = positions[id];
+	p.direction = directions[id];
+	p.direction /= norm(p.direction);
+	p.polarization = polarizations[id];
+	p.polarization /= norm(p.polarization);
+	p.wavelength = wavelengths[id];
+	p.time = times[id];
+	p.last_hit_triangle = last_hit_triangles[id];
+	p.history = histories[id];
+
+	if (p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB))
 		return;
 
-	curandState rng = rng_states[id];
-	float3 position = positions[id];
-	float3 direction = directions[id];
-	direction /= norm(direction);
-	float3 polarization = polarizations[id];
-	polarization /= norm(polarization);
-	float wavelength = wavelengths[id];
-	float time = times[id];
-	int last_hit_triangle = last_hit_triangles[id];
+	State s;
 
 	int steps = 0;
 	while (steps < max_steps)
 	{
 		steps++;
 
-		float distance;
+		int command;
 
-	        last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
+		command = fill_state(s, p);
 
-		if (last_hit_triangle == -1)
-		{
-			state = NO_HIT;
+		if (command == BREAK)
 			break;
-		}
-
-		uint4 triangle_data = triangles[last_hit_triangle];
-
-		float3 v0 = vertices[triangle_data.x];
-		float3 v1 = vertices[triangle_data.y];
-		float3 v2 = vertices[triangle_data.z];
-
-		int material_in_index = convert(0xFF & (triangle_data.w >> 24));
-		int material_out_index = convert(0xFF & (triangle_data.w >> 16));
-		int surface_index = convert(0xFF & (triangle_data.w >> 8));
-
-		float3 surface_normal = cross(v1-v0,v2-v1);
-		surface_normal /= norm(surface_normal);
-		
-		Material material1, material2;
-		if (dot(surface_normal,-direction) > 0.0f)
-		{
-			// outside to inside
-			material1 = materials[material_out_index];
-			material2 = materials[material_in_index];
-		}
-		else
-		{
-			// inside to outside
-			material1 = materials[material_in_index];
-			material2 = materials[material_out_index];
-			surface_normal = -surface_normal;
-		}
-
-		float refractive_index1 = interp_property(wavelength, material1.refractive_index);
-		float refractive_index2 = interp_property(wavelength, material2.refractive_index);
-
-		float absorption_length = interp_property(wavelength, material1.absorption_length);
-		float scattering_length = interp_property(wavelength, material1.scattering_length);
-
-		float absorption_distance = -absorption_length*logf(curand_uniform(&rng));
-		float scattering_distance = -scattering_length*logf(curand_uniform(&rng));
-
-		if (absorption_distance <= scattering_distance)
-		{
-			if (absorption_distance <= distance)
-			{
-				time += absorption_distance/(SPEED_OF_LIGHT/refractive_index1);
-				position += absorption_distance*direction;
-				state = BULK_ABSORB;
-
-				last_hit_triangle = -1;
 
-				break;
-			} // photon is absorbed in material1
-		}
-		else
-		{
-			if (scattering_distance <= distance)
-			{
-				time += scattering_distance/(SPEED_OF_LIGHT/refractive_index1);
-				position += scattering_distance*direction;
-
-				float theta,y;
-				while (true)
-				{
-					y = curand_uniform(&rng);
-					theta = uniform(&rng, 0, 2*PI);
-
-					if (y < powf(cosf(theta),2))
-						break;
-				}
-				
-				float phi = uniform(&rng, 0, 2*PI);
-
-				float3 b = cross(polarization, direction);
-				float3 c = polarization;
-
-				direction = rotate(direction, theta, b);
-				direction = rotate(direction, phi, c);
-				polarization = rotate(polarization, theta, b);
-				polarization = rotate(polarization, phi, c);
-
-				last_hit_triangle = -1;
-
-				continue;
-			} // photon is scattered in material1
-
-		} // if scattering_distance < absorption_distance
-
-		position += distance*direction;
-		time += distance/(SPEED_OF_LIGHT/refractive_index1);
+		command = propagate_to_boundary(p, s);
 
-		// p is normal to the plane of incidence
-		float3 p = cross(direction, surface_normal);
-		p /= norm(p);
-
-		float normal_coefficient = dot(polarization, p);
-		float normal_probability = normal_coefficient*normal_coefficient;
+		if (command == BREAK)
+			break;
 
-		float incident_angle = acosf(dot(surface_normal,-direction));
+		if (command == CONTINUE)
+			continue;
 
-		if (surface_index != -1)
+		if (s.surface_index != -1)
 		{
-			Surface surface = surfaces[surface_index];
-
-			float detect = interp_property(wavelength, surface.detect);
-			float absorb = interp_property(wavelength, surface.absorb);
-			float reflect_diffuse = interp_property(wavelength, surface.reflect_diffuse);
-			float reflect_specular = interp_property(wavelength, surface.reflect_specular);
-
-			// since the surface properties are interpolated linearly,
-			// we are guaranteed that they still sum to one.
-
-			float uniform_sample = curand_uniform(&rng);
+			command = propagate_at_surface(p, s);
 
-			if (uniform_sample < absorb)
-			{
-				// absorb
-				state = SURFACE_ABSORB;
+			if (command == BREAK)
 				break;
-			}
-			else if (uniform_sample < absorb + detect)
-			{
-				// detect
-				state = SURFACE_DETECT;
-				break;
-			}
-			else if (uniform_sample < absorb + detect + reflect_diffuse)
-			{
-				// diffusely reflect
-				direction = uniform_sphere(&rng);
-
-				if (dot(direction, surface_normal) < 0.0f)
-					direction = -direction;
-
-				// randomize polarization ?
-				polarization = cross(uniform_sphere(&rng), direction);
-				polarization /= norm(polarization);
 
+			if (command == CONTINUE)
 				continue;
-			}
-			else
-			{
-				// specularly reflect
-				direction = rotate(surface_normal, incident_angle, p);
-
-				// polarization ?
-
-				continue;
-			}
-
-		} // surface
-
-		float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2);
-
-		if (curand_uniform(&rng) < normal_probability)
-		{
-
-			float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
-			float reflection_probability = reflection_coefficient*reflection_coefficient;
-
-			if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
-			{
-				direction = rotate(surface_normal, incident_angle, p);
-			}
-			else
-			{
-				direction = rotate(surface_normal, PI-refracted_angle, p);
-			}
-
-			polarization = p;
-
-			continue;
-		}  // photon polarization normal to plane of incidence
-		else
-		{
-			float reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
-			float reflection_probability = reflection_coefficient*reflection_coefficient;
+		}
 
-			if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
-			{
-				direction = rotate(surface_normal, incident_angle, p);
-			}
-			else
-			{
-				direction = rotate(surface_normal, PI-refracted_angle, p);
-			}
+		propagate_at_boundary(p, s);
 
-			polarization = cross(p, direction);
-			polarization /= norm(polarization);
+	} // while (steps < max_steps)
 
-			continue;
-		} // photon polarization parallel to surface
-
-	} // while(nsteps < max_steps)
-
-	rng_states[id] = rng;
-	states[id] = state;
-	positions[id] = position;
-	directions[id] = direction;
-	polarizations[id] = polarization;
-	wavelengths[id] = wavelength;
-	times[id] = time;
-	last_hit_triangles[id] = last_hit_triangle;
+	rng_states[id] = p.rng;
+	positions[id] = p.position;
+	directions[id] = p.direction;
+	polarizations[id] = p.polarization;
+	wavelengths[id] = p.wavelength;
+	times[id] = p.time;
+	histories[id] = p.history;
+	last_hit_triangles[id] = p.last_hit_triangle;
 
 } // propagate
 
diff --git a/src/mesh.h b/src/mesh.h
new file mode 100644
index 0000000..6f678c0
--- /dev/null
+++ b/src/mesh.h
@@ -0,0 +1,146 @@
+#ifndef __MESH_H__
+#define __MESH_H__
+
+#include "intersect.h"
+
+#define STACK_SIZE 500
+
+/* flattened triangle mesh */
+__device__ float3 *g_vertices;
+__device__ uint4 *g_triangles;
+__device__ int g_start_node;
+__device__ int g_first_node;
+
+/* 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<unsigned int, 1, cudaReadModeElementType> node_map;
+texture<unsigned int, 1, cudaReadModeElementType> node_length;
+
+__device__ float3 make_float3(const float4 &a)
+{
+	return make_float3(a.x, a.y, a.z);
+}
+
+__device__ int convert(int c)
+{
+	if (c & 0x80)
+		return (0xFFFFFF00 | c);
+	else
+		return c;
+}
+
+/* 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, float &min_distance, int last_hit_triangle = -1)
+{
+	int triangle_index = -1;
+
+	float distance;
+
+	if (!intersect_node(origin, direction, g_start_node))
+		return -1;
+
+	int stack[STACK_SIZE];
+
+	int *head = &stack[0];
+	int *node = &stack[1];
+	int *tail = &stack[STACK_SIZE-1];
+	*node = g_start_node;
+
+	int i;
+
+	do
+	{
+		int first_child = tex1Dfetch(node_map, *node);
+		int child_count = tex1Dfetch(node_length, *node);
+
+		while (*node >= g_first_node && child_count == 1)
+		{
+			*node = first_child;
+			first_child = tex1Dfetch(node_map, *node);
+			child_count = tex1Dfetch(node_length, *node);
+		}
+		
+		if (*node >= g_first_node)
+		{
+			for (i=0; i < child_count; i++)
+			{
+				if (intersect_node(origin, direction, first_child+i))
+				{
+					*node = first_child+i;
+					node++;
+				}
+			}
+
+			node--;
+		}
+		else // node is a leaf
+		{
+			for (i=0; i < child_count; i++)
+			{
+				if (last_hit_triangle == first_child+i)
+					continue;
+
+				uint4 triangle_data = g_triangles[first_child+i];
+
+				float3 v0 = g_vertices[triangle_data.x];
+				float3 v1 = g_vertices[triangle_data.y];
+				float3 v2 = g_vertices[triangle_data.z];
+
+				if (intersect_triangle(origin, direction, v0, v1, v2, distance))
+				{
+					if (triangle_index == -1)
+					{
+						triangle_index = first_child + i;
+						min_distance = distance;
+						continue;
+					}
+
+					if (distance < min_distance)
+					{
+						triangle_index = first_child + i;
+						min_distance = distance;
+					}
+				}
+			} // triangle loop
+
+			node--;
+
+		} // node is a leaf
+
+	} // while loop
+	while (node != head);
+
+	return triangle_index;
+}
+
+extern "C"
+{
+
+__global__ void set_global_mesh_variables(uint4 *triangle_ptr, float3 *vertex_ptr, int start_node, int first_node)
+{
+	g_triangles = triangle_ptr;
+	g_vertices = vertex_ptr;
+	g_start_node = start_node;
+	g_first_node = first_node;
+}
+
+} // extern "c"
+
+#endif
diff --git a/src/photon.h b/src/photon.h
new file mode 100644
index 0000000..da866bc
--- /dev/null
+++ b/src/photon.h
@@ -0,0 +1,273 @@
+#ifndef __PHOTON_H__
+#define __PHOTON_H__
+
+#include "linalg.h"
+#include "materials.h"
+#include "rotate.h"
+#include "random.h"
+#include "physical_constants.h"
+#include "mesh.h"
+
+struct Photon
+{
+	float3 position;
+	float3 direction;
+	float3 polarization;
+	float wavelength;
+	float time;
+
+	//bool alive;
+	//unsigned int last_process;
+	unsigned int history;
+
+	int last_hit_triangle;
+
+	//int photon_state;
+
+	curandState rng;
+};
+
+struct State
+{
+	float3 surface_normal;
+	//Material material1, material2;
+	//Surface surface;
+
+	float refractive_index1, refractive_index2;
+	float absorption_length;
+	float scattering_length;
+	//float absorption_distance;
+	//float scattering_distance;
+
+	int surface_index;
+
+	//int last_hit_triangle;
+	float distance_to_boundary;
+};
+
+enum
+{
+	//DEBUG = -2,
+	//INIT = -1,
+	NO_HIT           = 0x1 << 0,
+	BULK_ABSORB      = 0x1 << 1,
+	SURFACE_DETECT   = 0x1 << 2,
+	SURFACE_ABSORB   = 0x1 << 3,
+	RAYLEIGH_SCATTER = 0x1 << 4,
+	REFLECT_DIFFUSE  = 0x1 << 5,
+	REFLECT_SPECULAR = 0x1 << 6
+}; // processes
+
+enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary
+
+__device__ int fill_state(State &s, Photon &p)
+{
+	p.last_hit_triangle = intersect_mesh(p.position, p.direction, s.distance_to_boundary, p.last_hit_triangle);
+
+	if (p.last_hit_triangle == -1)
+	{
+		p.history |= NO_HIT;
+		return BREAK;
+	}
+
+	uint4 triangle_data = g_triangles[p.last_hit_triangle];
+
+	float3 v0 = g_vertices[triangle_data.x];
+	float3 v1 = g_vertices[triangle_data.y];
+	float3 v2 = g_vertices[triangle_data.z];
+
+	int inner_material_index = convert(0xFF & (triangle_data.w >> 24));
+	int outer_material_index = convert(0xFF & (triangle_data.w >> 16));
+	s.surface_index = convert(0xFF & (triangle_data.w >> 8));
+
+	s.surface_normal = cross(v1-v0, v2-v1);
+	s.surface_normal /= norm(s.surface_normal);
+
+	Material material1, material2;
+	if (dot(s.surface_normal,-p.direction) > 0.0f)
+	{
+		// outside to inside
+		material1 = materials[outer_material_index];
+		material2 = materials[inner_material_index];
+	}
+	else
+	{
+		// inside to outside
+		material1 = materials[inner_material_index];
+		material2 = materials[outer_material_index];
+		s.surface_normal = -s.surface_normal;
+	}
+
+	s.refractive_index1 = interp_property(p.wavelength, material1.refractive_index);
+	s.refractive_index2 = interp_property(p.wavelength, material2.refractive_index);
+	s.absorption_length = interp_property(p.wavelength, material1.absorption_length);
+	s.scattering_length = interp_property(p.wavelength, material1.scattering_length);
+
+	return PASS;
+
+} // fill_state
+
+__device__ void rayleigh_scatter(Photon &p)
+{
+	float theta, y;
+
+	while (true)
+	{
+		y = curand_uniform(&p.rng);
+		theta = uniform(&p.rng, 0, 2*PI);
+
+		if (y < powf(cosf(theta),2))
+			break;
+	}
+
+	float phi = uniform(&p.rng, 0, 2*PI);
+
+	float3 b = cross(p.polarization, p.direction);
+	float3 c = p.polarization;
+
+	p.direction = rotate(p.direction, theta, b);
+	p.direction = rotate(p.direction, phi, c);
+
+	p.polarization = rotate(p.polarization, theta, b);
+	p.polarization = rotate(p.polarization, phi, c);
+
+} // scatter
+
+__device__ int propagate_to_boundary(Photon &p, State &s)
+{
+	float absorption_distance = -s.absorption_length*logf(curand_uniform(&p.rng));
+	float scattering_distance = -s.scattering_length*logf(curand_uniform(&p.rng));
+
+	if (absorption_distance <= scattering_distance)
+	{
+		if (absorption_distance <= s.distance_to_boundary)
+		{
+			p.time += absorption_distance/(SPEED_OF_LIGHT/s.refractive_index1);
+			p.position += absorption_distance*p.direction;
+			p.history |= BULK_ABSORB;
+
+			p.last_hit_triangle = -1;
+
+			return BREAK;
+		} // photon is absorbed in material1
+	}
+	else
+	{
+		if (scattering_distance <= s.distance_to_boundary)
+		{
+			p.time += scattering_distance/(SPEED_OF_LIGHT/s.refractive_index1);
+			p.position += scattering_distance*p.direction;
+
+			rayleigh_scatter(p);
+
+			p.history |= RAYLEIGH_SCATTER;
+
+			p.last_hit_triangle = -1;
+
+			return CONTINUE;
+		} // photon is scattered in material1
+	} // if scattering_distance < absorption_distance
+
+	p.position += s.distance_to_boundary*p.direction;
+	p.time += s.distance_to_boundary/(SPEED_OF_LIGHT/s.refractive_index1);
+
+	return PASS;
+
+} // propagate_to_boundary
+
+__device__ void propagate_at_boundary(Photon &p, State &s)
+{
+	float incident_angle = acosf(dot(s.surface_normal, -p.direction));
+	float refracted_angle = asinf(sinf(incident_angle)*s.refractive_index1/s.refractive_index2);
+
+	float3 incident_plane_normal = cross(p.direction, s.surface_normal);
+	incident_plane_normal /= norm(incident_plane_normal);
+
+	float normal_coefficient = dot(p.polarization, incident_plane_normal);
+	float normal_probability = normal_coefficient*normal_coefficient;
+
+	float reflection_coefficient;
+	if (curand_uniform(&p.rng) < normal_probability)
+	{
+		reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
+
+		p.polarization = incident_plane_normal;
+	} // photon polarization normal to plane of incidence
+	else
+	{
+		reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
+
+		p.polarization = cross(incident_plane_normal, p.direction);
+		p.polarization /= norm(p.polarization);
+	} // photon polarization parallel to plane of incidence
+
+	if ((curand_uniform(&p.rng) < reflection_coefficient*reflection_coefficient) || isnan(refracted_angle))
+	{
+		p.direction = rotate(s.surface_normal, incident_angle, incident_plane_normal);
+
+		p.history |= REFLECT_SPECULAR;
+	}
+	else
+	{
+		p.direction = rotate(s.surface_normal, PI-refracted_angle, incident_plane_normal);
+	}
+
+} // propagate_at_boundary
+
+__device__ int propagate_at_surface(Photon &p, State &s)
+{
+	Surface surface = surfaces[s.surface_index];
+
+	float detect = interp_property(p.wavelength, surface.detect);
+	float absorb = interp_property(p.wavelength, surface.absorb);
+	float reflect_diffuse = interp_property(p.wavelength, surface.reflect_diffuse);
+	float reflect_specular = interp_property(p.wavelength, surface.reflect_specular);
+
+	// since the surface properties are interpolated linearly, we are
+	// guaranteed that they still sum to 1.0.
+
+	float uniform_sample = curand_uniform(&p.rng);
+
+	if (uniform_sample < absorb)
+	{
+		p.history |= SURFACE_ABSORB;
+		return BREAK;
+	}
+	else if (uniform_sample < absorb + detect)
+	{
+		p.history |= SURFACE_DETECT;
+		return BREAK;
+	}
+	else if (uniform_sample < absorb + detect + reflect_diffuse)
+	{
+		// diffusely reflect
+		p.direction = uniform_sphere(&p.rng);
+
+		if (dot(p.direction, s.surface_normal) < 0.0f)
+			p.direction = -p.direction;
+
+		// randomize polarization?
+		p.polarization = cross(uniform_sphere(&p.rng), p.direction);
+		p.polarization /= norm(p.polarization);
+
+		p.history |= REFLECT_DIFFUSE;
+
+		return CONTINUE;
+	}
+	else
+	{
+		// specularly reflect
+		float incident_angle = acosf(dot(s.surface_normal, -p.direction));
+		float3 incident_plane_normal = cross(p.direction, s.surface_normal);
+		incident_plane_normal /= norm(incident_plane_normal);
+
+		p.direction = rotate(s.surface_normal, incident_angle, incident_plane_normal);
+
+		p.history |= REFLECT_SPECULAR;
+
+		return CONTINUE;
+	}
+
+} // propagate_at_surface
+
+#endif