when throwing photons from the light source out onto the scene, photons are now thrown randomly across each triangle instead of only at the center of each triangle. all of the rendering kernels have been rewritten so that they operate additively; for example, you may now throw photons from the light source onto the scene, render from the camera to the scene, then throw more photons and render again.

2 files changed, 100 insertions, 146 deletions

diff --git a/src/kernel.cu b/src/kernel.cu
index a36ff91..be0087d 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);
@@ -87,104 +85,76 @@ __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, float3 *positions, float3 *directions, float3 *rgb_lookup1, float3 *rgb_lookup2, unsigned long long seed, int runs, int max_steps)
+__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 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;
+	int kernel_id = blockIdx.x*blockDim.x + threadIdx.x;
+	int id = kernel_id + offset;
 
-	float distance;
-
-	int hit_triangle = intersect_mesh(p0.position, p0.direction, distance);
-
-	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;
 
-__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)
+} // update_xyz_lookup
+
+__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 +163,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 +280,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)