hybrid monte carlo render now distinguishes between the two different sides of each triangle. reduced the number of runs to average when propagating photons from each pixel in render.py from 5 to 1; the speed improvement outweighs any small improvement in the quality of the rendered image.

2 files changed, 113 insertions, 50 deletions

diff --git a/src/kernel.cu b/src/kernel.cu
index 73800a3..2cd9347 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -12,20 +12,14 @@
 #define BLUE_WAVELENGTH 465
 #define GREEN_WAVELENGTH 545
 
-__device__ void myAtomicAdd(float *addr, float data)
+__device__ void fAtomicAdd(float *addr, float data)
 {
 	while (data)
 		data = atomicExch(addr, data+atomicExch(addr, 0.0f));
 }
 
-__device__ int to_diffuse(Photon p, int max_steps)
+__device__ void to_diffuse(Photon &p, State &s, const int &max_steps)
 {
-	// 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)
 	{
@@ -51,7 +45,7 @@ __device__ int to_diffuse(Photon p, int max_steps)
 			command = propagate_at_surface(p, s);
 
 			if (p.history & REFLECT_DIFFUSE)
-				return p.last_hit_triangle;
+				break;
 
 			if (command == BREAK)
 				break;
@@ -64,7 +58,6 @@ __device__ int to_diffuse(Photon p, int max_steps)
 
 	} // while (steps < max_steps)
 
-	return -1;
 } // to_diffuse
 
 extern "C"
@@ -93,25 +86,25 @@ __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, float3 *rgb_lookup, int runs, int max_steps)
+__global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float3 *rgb_lookup1, float3 *rgb_lookup2, int runs, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	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;
+	Photon seed;
+	seed.rng = rng_states[id];
+	seed.position = positions[id];
+	seed.direction = directions[id];
+	seed.direction /= norm(seed.direction);
+	seed.polarization = uniform_sphere(&seed.rng);
+	seed.time = 0.0f;
+	seed.history = 0x0;
 
 	float distance;
 
-	int hit_triangle = intersect_mesh(p.position, p.direction, distance);
+	int hit_triangle = intersect_mesh(seed.position, seed.direction, distance);
 
 	if (hit_triangle != id)
 		return;
@@ -124,79 +117,143 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 	float3 v1 = g_vertices[triangle_data.y];
 	float3 v2 = g_vertices[triangle_data.z];
 
-	float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -p.direction);
+	float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -seed.direction);
 
 	if (cos_theta < 0.0f)
-		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -p.direction);
+		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -seed.direction);
+
+	Photon p;
+	State s;
 
 	for (int i=0; i < runs; i++)
 	{
+		p = seed;
 		p.wavelength = RED_WAVELENGTH;
 
-		hit_triangle = to_diffuse(p, max_steps);
+		to_diffuse(p, s, max_steps);
 
-		if (hit_triangle != -1)
-			myAtomicAdd(&rgb_lookup[hit_triangle].x, cos_theta);
+		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);
+			}
+		}
 
+		p = seed;
 		p.wavelength = BLUE_WAVELENGTH;
 
-		hit_triangle = to_diffuse(p, max_steps);
+		to_diffuse(p, s, max_steps);
 
-		if (hit_triangle != -1)
-			myAtomicAdd(&rgb_lookup[hit_triangle].y, cos_theta);
+		if (p.history & REFLECT_DIFFUSE)
+		{
+			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);
+			}
+		}
 
+		p = seed;
 		p.wavelength = GREEN_WAVELENGTH;
 
-		hit_triangle = to_diffuse(p, max_steps);
+		to_diffuse(p, s, max_steps);
 
-		if (hit_triangle != -1)
-			myAtomicAdd(&rgb_lookup[hit_triangle].z, cos_theta);
+		if (p.history & REFLECT_DIFFUSE)
+		{
+			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);
+			}
+		}
 	}
 
 } // build_rgb_lookup
 
-__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float3 *rgb_lookup, int runs, int *pixels, int max_steps)
+__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float3 *rgb_lookup1, float3 *rgb_lookup2, int runs, int *pixels, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (id >= nthreads)
 		return;
 
-	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;
+	Photon seed;
+	seed.rng = rng_states[id];
+	seed.position = positions[id];
+	seed.direction = directions[id];
+	seed.direction /= norm(seed.direction);
+	seed.polarization = uniform_sphere(&seed.rng);
+	seed.time = 0.0f;
+	seed.history = 0x0;
 
 	float3 rgb = make_float3(0.0, 0.0, 0.0);
 
-	int hit_triangle;
+	Photon p;
+	State s;
 
 	for (int i=0; i < runs; i++)
 	{
+		p = seed;
 		p.wavelength = RED_WAVELENGTH;
 
-		hit_triangle = to_diffuse(p, max_steps);
+		to_diffuse(p, s, max_steps);
 
-		if (hit_triangle != -1)
-			rgb.x += rgb_lookup[hit_triangle].x;
+		if (p.history & REFLECT_DIFFUSE)
+		{
+			if (s.inside_to_outside)
+			{
+				rgb.x += rgb_lookup1[p.last_hit_triangle].x;
+			}
+			else
+			{
+				rgb.x += rgb_lookup2[p.last_hit_triangle].x;
+			}
+		}
 
+		p = seed;
 		p.wavelength = BLUE_WAVELENGTH;
 
-		hit_triangle = to_diffuse(p, max_steps);
+		to_diffuse(p, s, max_steps);
 
-		if (hit_triangle != -1)
-			rgb.y += rgb_lookup[hit_triangle].y;
+		if (p.history & REFLECT_DIFFUSE)
+		{
+			if (s.inside_to_outside)
+			{
+				rgb.y += rgb_lookup1[p.last_hit_triangle].y;
+			}
+			else
+			{
+				rgb.y += rgb_lookup2[p.last_hit_triangle].y;
+			}
+		}
 
+		p = seed;
 		p.wavelength = GREEN_WAVELENGTH;
 
-		hit_triangle = to_diffuse(p, max_steps);
+		to_diffuse(p, s, max_steps);
 
-		if (hit_triangle != -1)
-			rgb.z += rgb_lookup[hit_triangle].z;
+		if (p.history & REFLECT_DIFFUSE)
+		{
+			if (s.inside_to_outside)
+			{
+				rgb.z += rgb_lookup1[p.last_hit_triangle].z;
+			}
+			else
+			{
+				rgb.z += rgb_lookup2[p.last_hit_triangle].z;
+			}
+		}
 	}
 
 	rgb /= runs;
diff --git a/src/photon.h b/src/photon.h
index 69ecd34..11bcfc1 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -25,6 +25,8 @@ struct Photon
 
 struct State
 {
+	bool inside_to_outside;
+
 	float3 surface_normal;
 
 	float refractive_index1, refractive_index2;
@@ -78,6 +80,8 @@ __device__ int fill_state(State &s, Photon &p)
 		// outside to inside
 		material1 = materials[outer_material_index];
 		material2 = materials[inner_material_index];
+
+		s.inside_to_outside = false;
 	}
 	else
 	{
@@ -85,6 +89,8 @@ __device__ int fill_state(State &s, Photon &p)
 		material1 = materials[inner_material_index];
 		material2 = materials[outer_material_index];
 		s.surface_normal = -s.surface_normal;
+
+		s.inside_to_outside = true;
 	}
 
 	s.refractive_index1 = interp_property(p.wavelength, material1.refractive_index);