pull random number generator states out of the photon struct; this allows you to copy photon information within a thread but still keep a single random number generator throughout the thread.

2 files changed, 46 insertions, 39 deletions

diff --git a/src/kernel.cu b/src/kernel.cu
index 2cd9347..711ae80 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -18,7 +18,7 @@ __device__ void fAtomicAdd(float *addr, float data)
 		data = atomicExch(addr, data+atomicExch(addr, 0.0f));
 }
 
-__device__ void to_diffuse(Photon &p, State &s, const int &max_steps)
+__device__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max_steps)
 {
 	int steps = 0;
 	while (steps < max_steps)
@@ -27,12 +27,12 @@ __device__ void to_diffuse(Photon &p, State &s, const int &max_steps)
 
 		int command;
 
-		command = fill_state(s, p);
+		command = fill_state(s, p, rng);
 
 		if (command == BREAK)
 			break;
 
-		command = propagate_to_boundary(p, s);
+		command = propagate_to_boundary(p, s, rng);
 
 		if (command == BREAK)
 			break;
@@ -42,7 +42,7 @@ __device__ void to_diffuse(Photon &p, State &s, const int &max_steps)
 
 		if (s.surface_index != -1)
 		{
-			command = propagate_at_surface(p, s);
+			command = propagate_at_surface(p, s, rng);
 
 			if (p.history & REFLECT_DIFFUSE)
 				break;
@@ -54,7 +54,7 @@ __device__ void to_diffuse(Photon &p, State &s, const int &max_steps)
 				continue;
 		}
 
-		propagate_at_boundary(p, s);
+		propagate_at_boundary(p, s, rng);
 
 	} // while (steps < max_steps)
 
@@ -93,12 +93,13 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 	if (id >= nthreads)
 		return;
 
+	curandState rng = rng_states[id];
+
 	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.polarization = uniform_sphere(&rng);
 	seed.time = 0.0f;
 	seed.history = 0x0;
 
@@ -130,7 +131,7 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 		p = seed;
 		p.wavelength = RED_WAVELENGTH;
 
-		to_diffuse(p, s, max_steps);
+		to_diffuse(p, s, rng, max_steps);
 
 		if (p.history & REFLECT_DIFFUSE)
 		{
@@ -147,7 +148,7 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 		p = seed;
 		p.wavelength = BLUE_WAVELENGTH;
 
-		to_diffuse(p, s, max_steps);
+		to_diffuse(p, s, rng, max_steps);
 
 		if (p.history & REFLECT_DIFFUSE)
 		{
@@ -164,7 +165,7 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 		p = seed;
 		p.wavelength = GREEN_WAVELENGTH;
 
-		to_diffuse(p, s, max_steps);
+		to_diffuse(p, s, rng, max_steps);
 
 		if (p.history & REFLECT_DIFFUSE)
 		{
@@ -179,6 +180,8 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 		}
 	}
 
+	rng_states[id] = rng;
+
 } // build_rgb_lookup
 
 __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)
@@ -188,12 +191,13 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 	if (id >= nthreads)
 		return;
 
+	curandState rng = rng_states[id];
+
 	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.polarization = uniform_sphere(&rng);
 	seed.time = 0.0f;
 	seed.history = 0x0;
 
@@ -207,7 +211,7 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 		p = seed;
 		p.wavelength = RED_WAVELENGTH;
 
-		to_diffuse(p, s, max_steps);
+		to_diffuse(p, s, rng, max_steps);
 
 		if (p.history & REFLECT_DIFFUSE)
 		{
@@ -224,7 +228,7 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 		p = seed;
 		p.wavelength = BLUE_WAVELENGTH;
 
-		to_diffuse(p, s, max_steps);
+		to_diffuse(p, s, rng, max_steps);
 
 		if (p.history & REFLECT_DIFFUSE)
 		{
@@ -241,7 +245,7 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 		p = seed;
 		p.wavelength = GREEN_WAVELENGTH;
 
-		to_diffuse(p, s, max_steps);
+		to_diffuse(p, s, rng, max_steps);
 
 		if (p.history & REFLECT_DIFFUSE)
 		{
@@ -264,6 +268,8 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 
 	pixels[id] = r << 16 | g << 8 | b;
 
+	rng_states[id] = rng;
+
 } // render
 
 /* Trace the rays starting at `positions` traveling in the direction
@@ -311,8 +317,9 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 	if (id >= nthreads)
 		return;
 
+	curandState rng = rng_states[id];
+
 	Photon p;
-	p.rng = rng_states[id];
 	p.position = positions[id];
 	p.direction = directions[id];
 	p.direction /= norm(p.direction);
@@ -335,12 +342,12 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 
 		int command;
 
-		command = fill_state(s, p);
+		command = fill_state(s, p, rng);
 
 		if (command == BREAK)
 			break;
 
-		command = propagate_to_boundary(p, s);
+		command = propagate_to_boundary(p, s, rng);
 
 		if (command == BREAK)
 			break;
@@ -350,7 +357,7 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 
 		if (s.surface_index != -1)
 		{
-			command = propagate_at_surface(p, s);
+			command = propagate_at_surface(p, s, rng);
 
 			if (command == BREAK)
 				break;
@@ -359,11 +366,11 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 				continue;
 		}
 
-		propagate_at_boundary(p, s);
+		propagate_at_boundary(p, s, rng);
 
 	} // while (steps < max_steps)
 
-	rng_states[id] = p.rng;
+	rng_states[id] = rng;
 	positions[id] = p.position;
 	directions[id] = p.direction;
 	polarizations[id] = p.polarization;
diff --git a/src/photon.h b/src/photon.h
index 11bcfc1..ad4c26c 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -20,7 +20,7 @@ struct Photon
 
 	int last_hit_triangle;
 
-	curandState rng;
+	//curandState rng;
 };
 
 struct State
@@ -51,7 +51,7 @@ enum
 
 enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary
 
-__device__ int fill_state(State &s, Photon &p)
+__device__ int fill_state(State &s, Photon &p, curandState &rng)
 {
 	p.last_hit_triangle = intersect_mesh(p.position, p.direction, s.distance_to_boundary, p.last_hit_triangle);
 
@@ -102,20 +102,20 @@ __device__ int fill_state(State &s, Photon &p)
 
 } // fill_state
 
-__device__ void rayleigh_scatter(Photon &p)
+__device__ void rayleigh_scatter(Photon &p, curandState &rng)
 {
 	float theta, y;
 
 	while (true)
 	{
-		y = curand_uniform(&p.rng);
-		theta = uniform(&p.rng, 0, 2*PI);
+		y = curand_uniform(&rng);
+		theta = uniform(&rng, 0, 2*PI);
 
 		if (y < powf(cosf(theta),2))
 			break;
 	}
 
-	float phi = uniform(&p.rng, 0, 2*PI);
+	float phi = uniform(&rng, 0, 2*PI);
 
 	float3 b = cross(p.polarization, p.direction);
 	float3 c = p.polarization;
@@ -128,10 +128,10 @@ __device__ void rayleigh_scatter(Photon &p)
 
 } // scatter
 
-__device__ int propagate_to_boundary(Photon &p, State &s)
+__device__ int propagate_to_boundary(Photon &p, State &s, curandState &rng)
 {
-	float absorption_distance = -s.absorption_length*logf(curand_uniform(&p.rng));
-	float scattering_distance = -s.scattering_length*logf(curand_uniform(&p.rng));
+	float absorption_distance = -s.absorption_length*logf(curand_uniform(&rng));
+	float scattering_distance = -s.scattering_length*logf(curand_uniform(&rng));
 
 	if (absorption_distance <= scattering_distance)
 	{
@@ -153,7 +153,7 @@ __device__ int propagate_to_boundary(Photon &p, State &s)
 			p.time += scattering_distance/(SPEED_OF_LIGHT/s.refractive_index1);
 			p.position += scattering_distance*p.direction;
 
-			rayleigh_scatter(p);
+			rayleigh_scatter(p, rng);
 
 			p.history |= RAYLEIGH_SCATTER;
 
@@ -170,7 +170,7 @@ __device__ int propagate_to_boundary(Photon &p, State &s)
 
 } // propagate_to_boundary
 
-__device__ void propagate_at_boundary(Photon &p, State &s)
+__device__ void propagate_at_boundary(Photon &p, State &s, curandState &rng)
 {
 	float incident_angle = acosf(dot(s.surface_normal, -p.direction));
 	float refracted_angle = asinf(sinf(incident_angle)*s.refractive_index1/s.refractive_index2);
@@ -182,11 +182,11 @@ __device__ void propagate_at_boundary(Photon &p, State &s)
 	float normal_probability = normal_coefficient*normal_coefficient;
 
 	float reflection_coefficient;
-	if (curand_uniform(&p.rng) < normal_probability)
+	if (curand_uniform(&rng) < normal_probability)
 	{
 		reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
 
-		if ((curand_uniform(&p.rng) < reflection_coefficient*reflection_coefficient) || isnan(refracted_angle))
+		if ((curand_uniform(&rng) < reflection_coefficient*reflection_coefficient) || isnan(refracted_angle))
 		{
 			p.direction = rotate(s.surface_normal, incident_angle, incident_plane_normal);
 			
@@ -203,7 +203,7 @@ __device__ void propagate_at_boundary(Photon &p, State &s)
 	{
 		reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
 
-		if ((curand_uniform(&p.rng) < reflection_coefficient*reflection_coefficient) || isnan(refracted_angle))
+		if ((curand_uniform(&rng) < reflection_coefficient*reflection_coefficient) || isnan(refracted_angle))
 		{
 			p.direction = rotate(s.surface_normal, incident_angle, incident_plane_normal);
 			
@@ -220,7 +220,7 @@ __device__ void propagate_at_boundary(Photon &p, State &s)
 
 } // propagate_at_boundary
 
-__device__ int propagate_at_surface(Photon &p, State &s)
+__device__ int propagate_at_surface(Photon &p, State &s, curandState &rng)
 {
 	Surface surface = surfaces[s.surface_index];
 
@@ -232,7 +232,7 @@ __device__ int propagate_at_surface(Photon &p, State &s)
 	// since the surface properties are interpolated linearly, we are
 	// guaranteed that they still sum to 1.0.
 
-	float uniform_sample = curand_uniform(&p.rng);
+	float uniform_sample = curand_uniform(&rng);
 
 	if (uniform_sample < absorb)
 	{
@@ -247,13 +247,13 @@ __device__ int propagate_at_surface(Photon &p, State &s)
 	else if (uniform_sample < absorb + detect + reflect_diffuse)
 	{
 		// diffusely reflect
-		p.direction = uniform_sphere(&p.rng);
+		p.direction = uniform_sphere(&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 = cross(uniform_sphere(&rng), p.direction);
 		p.polarization /= norm(p.polarization);
 
 		p.history |= REFLECT_DIFFUSE;