added a hybrid monte carlo ray tracing rendering algorithm

2 files changed, 440 insertions, 26 deletions

diff --git a/src/kernel.cu b/src/kernel.cu
index f257c6c..aae6e95 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -14,10 +14,12 @@
 
 enum
 {
+	DIFFUSE_HIT = -3,
 	DEBUG = -2,
 	INIT = -1,
 	NO_HIT,
 	BULK_ABSORB,
+	SURFACE_DETECT,
 	SURFACE_ABSORB
 };
 
@@ -144,28 +146,235 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 	return triangle_index;
 }
 
-extern "C"
+__device__ void myAtomicAdd(float *addr, float data)
 {
+	while (data)
+		data = atomicExch(addr, data+atomicExch(addr, 0.0f));
+}
+
+
 
-__global__ void fill_float(int nthreads, float *a, float value)
-{
-	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
-	if (id >= nthreads)
-		return;
 
-	a[id] = value;
-}
 
-__global__ void fill_float3(int nthreads, float3 *a, float3 value)
+
+__device__ int to_diffuse(curandState &rng, float3 position, float3 direction, float wavelength, float3 polarization, int start_node, int first_node, int max_steps)
 {
-	int id = blockIdx.x*blockDim.x + threadIdx.x;
+	int last_hit_triangle = -1;
+
+	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 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)
+			{
+				return -1;
+			} // 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);
+
+		// 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;
+
+		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 ?
+
+				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
+
+	} // while(nsteps < max_steps)
+
+	return -1;
+
+} // to_diffuse
+
+
 
-	if (id >= nthreads)
-		return;
 
-	a[id] = value;
-}
+
+
+
+
+
+
+
+
+
+
+
+
+
+extern "C"
+{
 
 __global__ void set_pointer(uint4 *triangle_ptr, float3 *vertex_ptr)
 {
@@ -196,6 +405,197 @@ __global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
 	points[id] = rotate(points[id], phi, axis);
 }
 
+#define RED_WAVELENGTH 685
+#define BLUE_WAVELENGTH 465
+#define GREEN_WAVELENGTH 545
+
+
+__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)
+{
+	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);
+
+	float distance;
+
+	int hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance);
+
+	if (hit_triangle != id)
+		return;
+
+	// note triangles from built geometry mesh are in order
+
+	uint4 triangle_data = triangles[hit_triangle];
+
+	float3 v0 = vertices[triangle_data.x];
+	float3 v1 = vertices[triangle_data.y];
+	float3 v2 = vertices[triangle_data.z];
+
+	float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -direction);
+
+	if (cos_theta < 0.0f)
+		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -direction);
+
+	for (int i=0; i < runs; i++)
+	{
+		hit_triangle = to_diffuse(rng, position, direction, RED_WAVELENGTH, polarization, start_node, first_node, 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);
+
+		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);
+
+		if (hit_triangle != -1)
+			myAtomicAdd(&rgb_lookup[hit_triangle].z, cos_theta);
+	}
+
+} // 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)
+{
+	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);
+
+	float3 rgb = make_float3(0.0, 0.0, 0.0);
+
+	int hit_triangle;
+
+	for (int i=0; i < runs; i++)
+	{
+		hit_triangle = to_diffuse(rng, position, direction, RED_WAVELENGTH, polarization, start_node, first_node, 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);
+
+		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);
+
+		if (hit_triangle != -1)
+			rgb.z += rgb_lookup[hit_triangle].z;
+	}
+
+	rgb /= runs;
+
+	unsigned int r = floorf(rgb.x*255);
+	unsigned int g = floorf(rgb.y*255);
+	unsigned int b = floorf(rgb.z*255);
+
+	pixels[id] = r << 16 | g << 8 | b;
+
+} // render
+
+#if 0
+
+__global__ void get_triangle(int nthreads, float3 *positions, float3 *directions, int start_node, int first_node, int *triangles)
+{
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (id >= nthreads)
+		return;
+
+	float3 position = positions[id];
+	float3 direction = directions[id];
+	direction /= norm(direction);
+
+	float distance;
+
+	triangles[id] = intersect_mesh(position, direction, start_node, first_node, distance);
+}
+
+__global__ void get_cos_theta(int nthreads, float3 *positions, float3 *directions, int start_node, int first_node, int *triangle, float *cos_thetas)
+{
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (id >= nthreads)
+		return;
+
+	float3 position = positions[id];
+	float3 direction = directions[id];
+	direction /= norm(direction);
+
+	uint4 triangle_data = triangles[triangle[id]];
+
+	float3 v0 = vertices[triangle_data.x];
+	float3 v1 = vertices[triangle_data.y];
+	float3 v2 = vertices[triangle_data.z];
+
+	float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -direction);
+
+	if (cos_theta < 0.0f)
+		cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -direction);
+
+	cos_thetas[id] = cos_theta;
+}
+
+#endif
+
+#if 0
+
+__global__ void to_diffuse(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, int *last_hit_triangles, int start_node, int first_node, int max_steps)
+{
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (id >= nthreads)
+		return;
+
+	curandState rng = rng_states[id];
+	float3 poisiton = positions[id];
+	float3 direction = directions[id];
+	direction /= norm(direction);
+	float3 polarization = polarizations[id];
+	polarization /= norm(polarization);
+	float wavelength = wavelengths[id];
+	float last_hit_triangle = last_hit_triangles[id];
+
+	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)
+			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);
+
+
+#endif
+
 /* 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
@@ -370,22 +770,29 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 		{
 			Surface surface = surfaces[surface_index];
 
-			float absorption = interp_property(wavelength, surface.absorption);
-			float reflection_diffuse = interp_property(wavelength, surface.reflection_diffuse);
-			float reflection_specular = interp_property(wavelength, surface.reflection_specular);
+			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 < absorption)
+			if (uniform_sample < absorb)
 			{
 				// absorb
 				state = SURFACE_ABSORB;
 				break;
 			}
-			else if (uniform_sample < absorption + reflection_diffuse)
+			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);
@@ -465,4 +872,9 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 
 } // propagate
 
+#if 0
+
+
+#endif
+
 } // extern "c"
diff --git a/src/materials.h b/src/materials.h
index 6115396..2355d24 100644
--- a/src/materials.h
+++ b/src/materials.h
@@ -15,9 +15,10 @@ struct Material
 
 struct Surface
 {
-	float *absorption;
-	float *reflection_diffuse;
-	float *reflection_specular;
+	float *detect;
+	float *absorb;
+	float *reflect_diffuse;
+	float *reflect_specular;
 };
 
 __device__ Material materials[20];
@@ -54,11 +55,12 @@ __global__ void set_material(int material_index, float *refractive_index, float
 	materials[material_index].scattering_length = scattering_length;
 }
 
-__global__ void set_surface(int surface_index, float *absorption, float *reflection_diffuse, float *reflection_specular)
+__global__ void set_surface(int surface_index, float *detect, float *absorb, float *reflect_diffuse, float *reflect_specular)
 {
-	surfaces[surface_index].absorption = absorption;
-	surfaces[surface_index].reflection_diffuse = reflection_diffuse;
-	surfaces[surface_index].reflection_specular = reflection_specular;
+	surfaces[surface_index].detect = detect;
+	surfaces[surface_index].absorb = absorb;
+	surfaces[surface_index].reflect_diffuse = reflect_diffuse;
+	surfaces[surface_index].reflect_specular = reflect_specular;
 }
 
 } // extern "c"