Merge with Tony

2 files changed, 112 insertions, 146 deletions

diff --git a/src/kernel.cu b/src/kernel.cu
index a36ff91..82efe88 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -9,14 +9,12 @@
 #include "photon.h"
 #include "alpha.h"
 
-#define RED_WAVELENGTH 685
-#define BLUE_WAVELENGTH 445
-#define GREEN_WAVELENGTH 545
-
 __device__ void fAtomicAdd(float *addr, float data)
 {
 	while (data)
+	{
 		data = atomicExch(addr, data+atomicExch(addr, 0.0f));
+	}
 }
 
 __device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max_steps)
@@ -28,9 +26,9 @@ __device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, c
 
 		int command;
 
-		command = fill_state(s, p, rng);
+		fill_state(s, p);
 
-		if (command == BREAK)
+		if (p.last_hit_triangle == -1)
 			break;
 
 		command = propagate_to_boundary(p, s, rng);
@@ -64,127 +62,111 @@ __device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, c
 extern "C"
 {
 
-/* Translate `points` by the vector `v` */
-__global__ void translate(int nthreads, float3 *points, float3 v)
+/* Translate the points `a` by the vector `v` */
+__global__ void translate(int nthreads, float3 *a, float3 v)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	points[id] += v;
+	a[id] += v;
 }
 
-/* Rotate `points` through an angle `phi` counter-clockwise about the
+/* Rotate the points `a` through an angle `phi` counter-clockwise about the
    axis `axis` (when looking towards +infinity). */
-__global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
+__global__ void rotate(int nthreads, float3 *a, float phi, float3 axis)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	points[id] = rotate(points[id], phi, axis);
+	a[id] = rotate(a[id], phi, axis);
 }
 
-  __global__ void build_rgb_lookup(int nthreads, float3 *positions, float3 *directions, float3 *rgb_lookup1, float3 *rgb_lookup2, unsigned long long seed, int runs, int max_steps)
+__global__ void rotate_around_point(int nthreads, float3 *a, float phi, float3 axis, float3 point)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	curandState rng;
-	curand_init(seed, id, 0, &rng);
-
-	Photon p0;
-	p0.position = positions[id];
-	p0.direction = directions[id];
-	p0.direction /= norm(p0.direction);
-	p0.polarization = uniform_sphere(&rng);
-	p0.time = 0.0f;
-	p0.history = 0x0;
-
-	float distance;
+	a[id] -= point;
+	a[id] = rotate(a[id], phi, axis);
+	a[id] += point;
+}
 
-	int hit_triangle = intersect_mesh(p0.position, p0.direction, distance);
+__global__ void update_xyz_lookup(int nthreads, int total_threads, int offset, float3 position, curandState *rng_states, float wavelength, float3 xyz, float3 *xyz_lookup1, float3 *xyz_lookup2, int max_steps)
+{
+	int kernel_id = blockIdx.x*blockDim.x + threadIdx.x;
+	int id = kernel_id + offset;
 
-	if (hit_triangle != id)
+	if (kernel_id >= nthreads || id >= total_threads)
 		return;
 
-	// note triangles from built geometry mesh are in order
+	curandState rng = rng_states[kernel_id];
 
-	uint4 triangle_data = g_triangles[hit_triangle];
+	uint4 triangle_data = g_triangles[id];
 
 	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)), -p0.direction);
+	float a = curand_uniform(&rng);
+	float b = uniform(&rng, 0.0f, (1.0f - a));
+	float c = 1.0f - a - b;
 
-	if (cos_theta < 0.0f)
-		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -p0.direction);
+	float3 direction = a*v0 + b*v1 + c*v2 - position;
+	direction /= norm(direction);
 
-	Photon p;
-	State s;
+	float distance;
+	int triangle_index = intersect_mesh(position, direction, distance);
 
-	for (int i=0; i < runs; i++)
+	if (triangle_index != id)
 	{
-		p = p0;
-		p.wavelength = RED_WAVELENGTH;
+		rng_states[kernel_id] = rng;
+		return;
+	}
 
-		to_diffuse(p, s, rng, max_steps);
+	float cos_theta = dot(normalize(cross(v1-v0,v2-v1)),-direction);
 
-		if (p.history & REFLECT_DIFFUSE)
-		{
-			if (s.inside_to_outside)
-			{
-				fAtomicAdd(&rgb_lookup1[p.last_hit_triangle].x, cos_theta);
-			}
-			else
-			{
-				fAtomicAdd(&rgb_lookup2[p.last_hit_triangle].x, cos_theta);
-			}
-		}
+	if (cos_theta < 0.0f)
+		cos_theta = dot(-normalize(cross(v1-v0,v2-v1)),-direction);
 
-		p = p0;
-		p.wavelength = GREEN_WAVELENGTH;
+	Photon p;
+	p.position = position;
+	p.direction = direction;
+	p.wavelength = wavelength;
+	p.polarization = uniform_sphere(&rng);
+	p.last_hit_triangle = -1;
+	p.time = 0;
+	p.history = 0;
 
-		to_diffuse(p, s, rng, max_steps);
+	State s;
+	to_diffuse(p, s, rng, max_steps);
 
-		if (p.history & REFLECT_DIFFUSE)
+	if (p.history & REFLECT_DIFFUSE)
+	{
+		if (s.inside_to_outside)
 		{
-			if (s.inside_to_outside)
-			{
-				fAtomicAdd(&rgb_lookup1[p.last_hit_triangle].y, cos_theta);
-			}
-			else
-			{
-				fAtomicAdd(&rgb_lookup2[p.last_hit_triangle].y, cos_theta);
-			}
+			fAtomicAdd(&xyz_lookup1[p.last_hit_triangle].x, cos_theta*xyz.x);
+			fAtomicAdd(&xyz_lookup1[p.last_hit_triangle].y, cos_theta*xyz.y);
+			fAtomicAdd(&xyz_lookup1[p.last_hit_triangle].z, cos_theta*xyz.z);
 		}
-
-		p = p0;
-		p.wavelength = BLUE_WAVELENGTH;
-
-		to_diffuse(p, s, rng, max_steps);
-
-		if (p.history & REFLECT_DIFFUSE)
+		else
 		{
-			if (s.inside_to_outside)
-			{
-				fAtomicAdd(&rgb_lookup1[p.last_hit_triangle].z, cos_theta);
-			}
-			else
-			{
-				fAtomicAdd(&rgb_lookup2[p.last_hit_triangle].z, cos_theta);
-			}
+			fAtomicAdd(&xyz_lookup2[p.last_hit_triangle].x, cos_theta*xyz.x);
+			fAtomicAdd(&xyz_lookup2[p.last_hit_triangle].y, cos_theta*xyz.y);
+			fAtomicAdd(&xyz_lookup2[p.last_hit_triangle].z, cos_theta*xyz.z);
 		}
 	}
 
-} // build_rgb_lookup
+	rng_states[kernel_id] = rng;
+
+} // update_xyz_lookup
 
-__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *red, float *green, float *blue, float3 *rgb_lookup1, float3 *rgb_lookup2, int max_steps)
+__global__ void update_xyz_image(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float wavelength, float3 xyz, float3 *xyz_lookup1, float3 *xyz_lookup2, float3 *image, int nlookup_calls, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
@@ -193,81 +175,69 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 
 	curandState rng = rng_states[id];
 
-	Photon p0;
-	p0.position = positions[id];
-	p0.direction = directions[id];
-	p0.direction /= norm(p0.direction);
-	p0.polarization = uniform_sphere(&rng);
-	p0.time = 0.0f;
-	p0.history = 0x0;
+	Photon p;
+	p.position = positions[id];
+	p.direction = directions[id];
+	p.direction /= norm(p.direction);
+	p.wavelength = wavelength;
+	p.polarization = uniform_sphere(&rng);
+	p.last_hit_triangle = -1;
+	p.time = 0;
+	p.history = 0;
 
 	State s;
-	Photon p = p0;
-	p.wavelength = RED_WAVELENGTH;
-		
 	to_diffuse(p, s, rng, max_steps);
 
 	if (p.history & REFLECT_DIFFUSE)
 	{
 		if (s.inside_to_outside)
 		{
-			red[id] = rgb_lookup1[p.last_hit_triangle].x;
-		}
-		else
-		{
-			red[id] = rgb_lookup2[p.last_hit_triangle].x;
-		}
-	}
-	else
-	{
-		red[id] = 0.0;
-	}
-	
-	p = p0;
-	p.wavelength = GREEN_WAVELENGTH;
-	
-	to_diffuse(p, s, rng, max_steps);
-	
-	if (p.history & REFLECT_DIFFUSE)
-	{
-		if (s.inside_to_outside)
-		{
-			green[id] = rgb_lookup1[p.last_hit_triangle].y;
+			image[id].x += xyz.x*xyz_lookup1[p.last_hit_triangle].x/nlookup_calls;
+			image[id].y += xyz.y*xyz_lookup1[p.last_hit_triangle].y/nlookup_calls;
+			image[id].z += xyz.z*xyz_lookup1[p.last_hit_triangle].z/nlookup_calls;
 		}
 		else
 		{
-			green[id] = rgb_lookup2[p.last_hit_triangle].y;
+			image[id].x += xyz.x*xyz_lookup2[p.last_hit_triangle].x/nlookup_calls;
+			image[id].y += xyz.y*xyz_lookup2[p.last_hit_triangle].y/nlookup_calls;
+			image[id].z += xyz.z*xyz_lookup2[p.last_hit_triangle].z/nlookup_calls;
 		}
 	}
-	else
-	{
-		green[id] = 0.0;
-	}
-	
-	p = p0;
-	p.wavelength = BLUE_WAVELENGTH;
-	
-	to_diffuse(p, s, rng, max_steps);
-	
-	if (p.history & REFLECT_DIFFUSE)
-	{
-		if (s.inside_to_outside)
-		{
-			blue[id] = rgb_lookup1[p.last_hit_triangle].z;
-		}
-		else
-		{
-			blue[id] = rgb_lookup2[p.last_hit_triangle].z;
-		}
-	}
-	else
-	{
-		blue[id] = 0.0;
-	}
-	
+
 	rng_states[id] = rng;
 
-} // render
+} // update_xyz_image
+
+__global__ void process_image(int nthreads, float3 *image, int *pixels, int nimages)
+{
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (id >= nthreads)
+		return;
+
+	float3 rgb = image[id]/nimages;
+
+	if (rgb.x < 0.0f)
+		rgb.x = 0.0f;
+	if (rgb.y < 0.0f)
+		rgb.y = 0.0f;
+	if (rgb.z < 0.0f)
+		rgb.z = 0.0f;
+
+	if (rgb.x > 1.0f)
+		rgb.x = 1.0f;
+	if (rgb.y > 1.0f)
+		rgb.y = 1.0f;
+	if (rgb.z > 1.0f)
+		rgb.z = 1.0f;
+
+	unsigned int r = floorf(rgb.x*255.0f);
+	unsigned int g = floorf(rgb.y*255.0f);
+	unsigned int b = floorf(rgb.z*255.0f);
+
+	pixels[id] = r << 16 | g << 8 | b;
+
+} // process_image
 
 /* Trace the rays starting at `positions` traveling in the direction
    `directions` to their intersection with the global mesh. If the ray
@@ -322,9 +292,9 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 
 		int command;
 
-		command = fill_state(s, p, rng);
+		fill_state(s, p);
 
-		if (command == BREAK)
+		if (p.last_hit_triangle == -1)
 			break;
 
 		command = propagate_to_boundary(p, s, rng);
diff --git a/src/photon.h b/src/photon.h
index ad4c26c..f471866 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -19,8 +19,6 @@ struct Photon
 	unsigned int history;
 
 	int last_hit_triangle;
-
-	//curandState rng;
 };
 
 struct State
@@ -51,14 +49,14 @@ enum
 
 enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary
 
-__device__ int fill_state(State &s, Photon &p, curandState &rng)
+__device__ void 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;
+		return;
 	}
 
 	uint4 triangle_data = g_triangles[p.last_hit_triangle];
@@ -98,8 +96,6 @@ __device__ int fill_state(State &s, Photon &p, curandState &rng)
 	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, curandState &rng)