added some fun models; added some untested code to implement absorption, scattering, reflection, and refraction

5 files changed, 318 insertions, 30 deletions

diff --git a/src/intersect.h b/src/intersect.h
index 92bcf0c..e6566f3 100644
--- a/src/intersect.h
+++ b/src/intersect.h
@@ -1,4 +1,8 @@
 //-*-c-*-
+
+#ifndef __INTERSECT_H__
+#define __INTERSECT_H__
+
 #include <math_constants.h>
 
 #include "linalg.h"
@@ -145,3 +149,5 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con
 
 	return true;
 }
+
+#endif
diff --git a/src/kernel.cu b/src/kernel.cu
index 796de54..b2e0903 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -7,18 +7,20 @@
 #include "rotate.h"
 #include "intersect.h"
 #include "materials.h"
+#include "physical_constants.h"
+#include "random.h"
 
 #define STACK_SIZE 500
+#define MAX_DEPTH 5
 
 /* flattened triangle mesh */
 texture<float4, 1, cudaReadModeElementType> vertices;
 __device__ uint4 *triangles;
 
 /* material/surface index lookup for each triangle */
-texture<int, 1, cudaReadModeElementType> material1_index;
-texture<int, 1, cudaReadModeElementType> material2_index;
-texture<int, 1, cudaReadModeElementType> surface1_index;
-texture<int, 1, cudaReadModeElementType> surface2_index;
+texture<int, 1, cudaReadModeElementType> material1_lookup;
+texture<int, 1, cudaReadModeElementType> material2_lookup;
+texture<int, 1, cudaReadModeElementType> surface_lookup;
 
 /* lower/upper bounds for the bounding box associated with each node/leaf */
 texture<float4, 1, cudaReadModeElementType> upper_bounds;
@@ -48,11 +50,10 @@ __device__ bool intersect_node(const float3 &origin, const float3 &direction, co
    direction `direction` and the global mesh texture. If the ray intersects
    the texture return the index of the triangle which the ray intersected,
    else return -1. */
-__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node)
+__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node, float &distance)
 {
 	int triangle_idx = -1;
 
-	float distance;
 	float min_distance;
 
 	if (!intersect_node(origin, direction, start_node))
@@ -133,20 +134,8 @@ __global__ void set_pointer(uint4 *triangle_ptr)
 /* Initialize random number states */
 __global__ void init_rng(unsigned long long seed, unsigned long long offset)
 {
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
-	curand_init(seed, idx, offset, rng_states+idx);
-}
-
-/* */
-__global__ void uniform_sphere(int nthreads, float3 *points)
-{
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
-
-	if (idx >= nthreads)
-		return;
-
-	//curandState rng = *(rng_states+idx);
-
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
+	curand_init(seed, id, offset, rng_states+id);
 }
 
 /* Translate `points` by the vector `v` */
@@ -188,7 +177,9 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int
 	float3 direction = *(directions+idx);
 	direction /= norm(direction);
 
-	int intersection_idx = intersect_mesh(origin, direction, start_node, first_node);
+	float distance;
+
+	int intersection_idx = intersect_mesh(origin, direction, start_node, first_node, distance);
 
 	if (intersection_idx == -1)
 	{
@@ -206,23 +197,244 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int
 	}
 } // ray_trace
 
-/* Propagate the photons starting at `origins` traveling in the direction
+/* Propagate the photons starting at `positions` traveling in the direction
    `directions` to their intersection with the global mesh. If the ray
    intersects the mesh set the hit_solid array value associated with the
    photon to the triangle index of the triangle the photon intersected, else
    set the hit_solid array value to -1. */
-__global__ void propagate(int nthreads, float3 *origins, float3 *directions, int start_node, int first_node, int *hit_triangles)
+__global__ void propagate(int nthreads, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, int *states, int *last_hit_triangles, int start_node, int first_node)
 {
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
-	if (idx >= nthreads)
+	//states[id] = interp_property(materials[0].refractive_index
+
+	//return;
+
+	if (id >= nthreads)
 		return;
 
-	float3 origin = *(origins+idx);
-	float3 direction = *(directions+idx);
+	float3 position = positions[id];
+	float3 direction = directions[id];
 	direction /= norm(direction);
+	float3 polarization = polarizations[id];
+	float wavelength = wavelengths[id];
+	float time = times[id];
+
+	int last_hit_triangle;
+
+	curandState rng = rng_states[id];
+
+	unsigned int depth=0;
+	while (true)
+	{
+		depth++;
+
+		if (depth > MAX_DEPTH)
+		{
+			states[id] = MAX_DEPTH_REACHED;
+			break;
+		}
+
+		float distance;
+		
+	        last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance);
+		
+		if (last_hit_triangle == -1)
+		{
+			states[id] = NO_HIT;
+			break;
+		}
+
+		uint4 triangle_data = triangles[last_hit_triangle];
+
+		float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x));
+		float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y));
+		float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z));
+		
+		int material_in_index = 0xFF & (triangle_data.w >> 24);
+		int material_out_index = 0xFF & (triangle_data.w >> 16);
+		int surface_index = 0xFF & (triangle_data.w >> 8);
+
+		if (material_in_index < 0 || material_out_index < 0 || surface_index < 0)
+		{
+			states[id] = -2;
+			break;
+		}
+		
+		float3 surface_normal = cross(v1-v0,v2-v0);
+		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 = -absorption_length*logf(curand_uniform(&rng));
+
+		if (absorption_distance <= scattering_distance)
+		{
+			if (absorption_distance <= distance)
+			{
+				time = absorption_distance/(SPEED_OF_LIGHT/refractive_index1);
+				position = position + absorption_distance*direction;
+				states[id] = BULK_ABSORB;
+				break;
+			} // photon is absorbed in material1
+		}
+		else
+		{
+			if (scattering_distance <= distance)
+			{
+				time = scattering_distance/(SPEED_OF_LIGHT/refractive_index1);
+				position = position + scattering_distance*direction;
+
+				float x,y;
+				while (true)
+				{
+					y = curand_uniform(&rng);
+					x = uniform(0, 2*PI);
+
+					if (y < powf(cosf(x),2))
+						break;
+				}
+				
+				float phi = uniform(0, 2*PI);
+
+				float3 old_direction = direction;
+				float3 old_polarization = polarization;
+
+				direction = rotate(old_direction, x, cross(old_polarization, old_direction));
+				polarization = rotate(old_polarization, x, cross(old_polarization, old_direction));
+
+				direction = rotate(direction, phi, old_polarization);
+				polarization = rotate(polarization, phi, old_polarization);
+
+				direction /= norm(direction);
+				polarization /= norm(polarization);
+
+				continue;
+			} // photon is scattered in material1
+		}
+		
+		if (surface_index != -1)
+		{
+			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 sum = absorption + reflection_diffuse + reflection_specular;
+			absorption /= sum;
+			reflection_diffuse /= sum;
+			reflection_specular /= sum;
+
+			float p = curand_uniform(&rng);
+
+			if (p < absorption)
+			{
+				// absorb
+				states[id] = SURFACE_ABSORB;
+				break;
+			}
+			else if (p < absorption + reflection_diffuse)
+			{
+				// diffusely reflect
+				while (true)
+				{
+					direction = uniform_sphere(&rng);
+
+					if (dot(direction, surface_normal) > 0.0f)
+					{
+						// randomize polarization
+						polarization = cross(uniform_sphere(&rng), direction);
+						break;
+					}
+				}
+
+				continue;
+			}
+			else
+			{
+				// specularly reflect
+				direction = rotate(direction, -PI/2.0f, cross(direction, surface_normal));
+
+				// polarization ?
+
+				continue;
+			}
+		} // surface
+		
+		// compute reflection/transmission probability
+
+		// p is normal to the plane of incidence
+		float3 p = cross(direction, surface_normal);
+
+		float normal_coefficient = dot(polarization, p);
+		float normal_probability = normal_coefficient*normal_coefficient;
+
+		float incident_angle = acosf(dot(surface_normal,-direction));
+		float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2);
+
+		if (curand_uniform(&rng) < normal_probability)
+		{
+			// photon polarization normal to plane of incidence
+			float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
+			float reflection_probability = reflection_coefficient*reflection_coefficient;
+
+			polarization = p;
+
+			if (curand_uniform(&rng) < reflection_probability)
+			{
+				// reflect off surface
+				direction = rotate(surface_normal, PI/2.0f, p);
+			}
+			else
+			{
+				// transmit
+				direction = rotate(surface_normal, PI-refracted_angle, p);
+			}
+			continue;
+		}
+		else
+		{
+			// photon polarization parallel to surface
+			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)
+			{
+				// reflect off surface
+				direction = rotate(surface_normal, PI/2.0f, p);
+				polarization = cross(p, direction);
+			}
+			else
+			{
+				// transmit
+				direction = rotate(surface_normal, PI-refracted_angle, p);
+				polarization = cross(p, direction);
+			}
+			continue;
+		}
 
-	*(hit_triangles+idx) = intersect_mesh(origin, direction, start_node, first_node);
+	} // while(true)
 
 } // propagate
 
diff --git a/src/materials.h b/src/materials.h
index 05aa121..47b7d22 100644
--- a/src/materials.h
+++ b/src/materials.h
@@ -1,3 +1,20 @@
+#ifndef __MATERIALS_H__
+#define __MATERIALS_H__
+
+__device__ float min_wavelength;
+__device__ float max_wavelength;
+__device__ float wavelength_step;
+__device__ unsigned int wavelength_size;
+
+enum
+{
+	INIT = -1,
+	MAX_DEPTH_REACHED,
+	NO_HIT,
+	BULK_ABSORB,
+	SURFACE_ABSORB
+};
+
 struct Material
 {
 	float *refractive_index;
@@ -12,12 +29,34 @@ struct Surface
 	float *reflection_specular;
 };
 
-__device__ struct Material materials[100];
-__device__ struct Surface surfaces[100];
+__device__ Material materials[100];
+__device__ Surface surfaces[100];
+
+__device__ float interp_property(const float &x, const float *fp)
+{
+	if (x < min_wavelength)
+		return fp[0];
+
+	if (x > max_wavelength)
+		return fp[wavelength_size-1];
+
+	unsigned int jl = (x-min_wavelength)/wavelength_step;
+	float xl = min_wavelength + jl*wavelength_step;
+
+	return xl + (x-xl)*(fp[jl+1]-fp[jl])/wavelength_step;
+}
 
 extern "C"
 {
 
+__global__ void set_wavelength_range(float _min_wavelength, float _max_wavelength, float _wavelength_step, unsigned int _wavelength_size)
+{
+	min_wavelength = _min_wavelength;
+	max_wavelength = _max_wavelength;
+	wavelength_step = _wavelength_step;
+	wavelength_size = _wavelength_size;
+}
+
 __global__ void set_material(int material_index, float *refractive_index, float *absorption_length, float *scattering_length)
 {
 	materials[material_index].refractive_index = refractive_index;
@@ -33,3 +72,5 @@ __global__ void set_surface(int surface_index, float *absorption, float *reflect
 }
 
 } // extern "c"
+
+#endif
diff --git a/src/physical_constants.h b/src/physical_constants.h
new file mode 100644
index 0000000..2ff87cd
--- /dev/null
+++ b/src/physical_constants.h
@@ -0,0 +1,7 @@
+#ifndef __PHYSICAL_CONSTANTS_H__
+#define __PHYSICAL_CONSTANTS_H__
+
+#define SPEED_OF_LIGHT 2.99792458e8
+#define PI 3.141592653589793
+
+#endif
diff --git a/src/random.h b/src/random.h
new file mode 100644
index 0000000..74329da
--- /dev/null
+++ b/src/random.h
@@ -0,0 +1,22 @@
+#ifndef __RANDOM_H__
+#define __RANDOM_H__
+
+#include <curand_kernel.h>
+
+#include "physical_constants.h"
+
+__device__ float uniform(curandState *rng, float a=0.0, float b=1.0)
+{
+	return a + curand_uniform(rng)*(b-a);
+}
+
+__device__ float3 uniform_sphere(curandState *rng)
+{
+	float theta = uniform(rng, 0, 2*PI);
+	float u = uniform(rng, -1, 1);
+	float c = sqrtf(1-u*u);
+
+	return make_float3(c*cosf(theta), c*sinf(theta), u);
+}
+
+#endif