1 files changed, 0 insertions, 389 deletions
diff --git a/src/kernel.cu b/src/kernel.cu
deleted file mode 100644
index 4418305..0000000
--- a/src/kernel.cu
+++ /dev/null
@@ -1,389 +0,0 @@
-//-*-c-*-
-#include <math_constants.h>
-#include <curand_kernel.h>
-
-#include "linalg.h"
-#include "matrix.h"
-#include "rotate.h"
-#include "mesh.h"
-#include "photon.h"
-#include "alpha.h"
-
-__device__ void fAtomicAdd(float *addr, float data)
-{
-	while (data)
-	{
-		data = atomicExch(addr, data+atomicExch(addr, 0.0f));
-	}
-}
-
-__device__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max_steps)
-{
-	int steps = 0;
-	while (steps < max_steps)
-	{
-		steps++;
-
-		int command;
-
-		fill_state(s, p);
-
-		if (p.last_hit_triangle == -1)
-			break;
-
-		command = propagate_to_boundary(p, s, rng);
-
-		if (command == BREAK)
-			break;
-
-		if (command == CONTINUE)
-			continue;
-
-		if (s.surface_index != -1)
-		{
-			command = propagate_at_surface(p, s, rng);
-
-			if (p.history & REFLECT_DIFFUSE)
-				break;
-
-			if (command == BREAK)
-				break;
-
-			if (command == CONTINUE)
-				continue;
-		}
-
-		propagate_at_boundary(p, s, rng);
-
-	} // while (steps < max_steps)
-
-} // to_diffuse
-
-extern "C"
-{
-
-__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 kernel_id = blockIdx.x*blockDim.x + threadIdx.x;
-	int id = kernel_id + offset;
-
-	if (kernel_id >= nthreads || id >= total_threads)
-		return;
-
-	curandState rng = rng_states[kernel_id];
-
-	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 a = curand_uniform(&rng);
-	float b = uniform(&rng, 0.0f, (1.0f - a));
-	float c = 1.0f - a - b;
-
-	float3 direction = a*v0 + b*v1 + c*v2 - position;
-	direction /= norm(direction);
-
-	float distance;
-	int triangle_index = intersect_mesh(position, direction, distance);
-
-	if (triangle_index != id)
-	{
-		rng_states[kernel_id] = rng;
-		return;
-	}
-
-	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);
-
-	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;
-
-	State s;
-	to_diffuse(p, s, rng, max_steps);
-
-	if (p.history & REFLECT_DIFFUSE)
-	{
-		if (s.inside_to_outside)
-		{
-			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);
-		}
-		else
-		{
-			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);
-		}
-	}
-
-	rng_states[kernel_id] = rng;
-
-} // 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;
-
-	if (id >= nthreads)
-		return;
-
-	curandState rng = rng_states[id];
-
-	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;
-	to_diffuse(p, s, rng, max_steps);
-
-	if (p.history & REFLECT_DIFFUSE)
-	{
-		if (s.inside_to_outside)
-		{
-			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
-		{
-			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;
-		}
-	}
-
-	rng_states[id] = rng;
-
-} // 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
-
-__global__ void distance_to_mesh(int nthreads, float3 *positions, float3 *directions, float *distances)
-{
-	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;
-
-	int triangle_index = intersect_mesh(position, direction, distance);
-
-	if (triangle_index == -1)
-		distances[id] = 1e9;
-	else
-		distances[id] = distance;
-
-} // distance_to_mesh
-
-__global__ void ray_trace(int nthreads, float3 *positions, float3 *directions, int *pixels)
-{
-	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;
-
-	int triangle_index = intersect_mesh(position, direction, distance);
-
-	if (triangle_index == -1)
-	{
-		pixels[id] = 0;
-	}
-	else
-	{
-		uint4 triangle_data = g_triangles[triangle_index];
-
-		float3 v0 = g_vertices[triangle_data.x];
-		float3 v1 = g_vertices[triangle_data.y];
-		float3 v2 = g_vertices[triangle_data.z];
-
-		pixels[id] = get_color(direction, v0, v1, v2, g_colors[triangle_index]);
-	}
-
-} // ray_trace
-
-/* 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 determined by the cosine of the angle between
-   the ray and the normal of the triangle it intersected, else set the pixel
-   to 0. */
-__global__ void ray_trace_alpha(int nthreads, float3 *positions, float3 *directions, int *pixels)
-{
-	int id = blockIdx.x*blockDim.x + threadIdx.x;
-
-	if (id >= nthreads)
-		return;
-
-	float3 position = positions[id];
-	float3 direction = directions[id];
-	direction /= norm(direction);
-
-	bool hit;
-	float distance;
-
-	pixels[id] = get_color_alpha(position, direction);
-
-} // ray_trace
-
-__global__ void swap(int *values, int nswap, int *offset_a, int *offset_b)
-{
-  int id = blockIdx.x*blockDim.x + threadIdx.x;
-
-  if (id < nswap) {
-    int a = offset_a[id];
-    int b = offset_b[id];
-
-    int tmp = values[a];
-    values[a] = values[b];
-    values[b] = tmp;
-  }
-}
-
-__global__ void propagate(int first_photon, int nthreads, 
-			  unsigned int *input_queue, 
-			  unsigned int *output_queue,
-			  curandState *rng_states, 
-			  float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, 
-			  unsigned int *histories, int *last_hit_triangles, int max_steps)
-{
-	int id = blockIdx.x*blockDim.x + threadIdx.x;
-
-	if (id >= nthreads)
-		return;
-
-	curandState rng = rng_states[id];
-
-	int photon_id = input_queue[first_photon + id];
-
-	Photon p;
-	p.position = positions[photon_id];
-	p.direction = directions[photon_id];
-	p.direction /= norm(p.direction);
-	p.polarization = polarizations[photon_id];
-	p.polarization /= norm(p.polarization);
-	p.wavelength = wavelengths[photon_id];
-	p.time = times[photon_id];
-	p.last_hit_triangle = last_hit_triangles[photon_id];
-	p.history = histories[photon_id];
-
-	if (p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB))
-		return;
-
-	State s;
-
-	int steps = 0;
-	while (steps < max_steps)
-	{
-		steps++;
-
-		int command;
-
-		// check for NaN and fail
-		if (isnan(p.direction.x*p.direction.y*p.direction.z*p.position.x*p.position.y*p.position.z)) 
-		{
-			p.history |= NO_HIT | NAN_ABORT;
-			break;
-		}
-
-		fill_state(s, p);
-
-		if (p.last_hit_triangle == -1)
-			break;
-
-		command = propagate_to_boundary(p, s, rng);
-
-		if (command == BREAK)
-			break;
-
-		if (command == CONTINUE)
-			continue;
-
-		if (s.surface_index != -1)
-		{
-			command = propagate_at_surface(p, s, rng);
-
-			if (command == BREAK)
-				break;
-
-			if (command == CONTINUE)
-				continue;
-		}
-
-		propagate_at_boundary(p, s, rng);
-
-	} // while (steps < max_steps)
-
-	rng_states[id] = rng;
-	positions[photon_id] = p.position;
-	directions[photon_id] = p.direction;
-	polarizations[photon_id] = p.polarization;
-	wavelengths[photon_id] = p.wavelength;
-	times[photon_id] = p.time;
-	histories[photon_id] = p.history;
-	last_hit_triangles[photon_id] = p.last_hit_triangle;
-	
-	// Not done, put photon in output queue
-	if ( (p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB)) == 0) {
-	  int out_idx = atomicAdd(output_queue, 1);
-	  output_queue[out_idx] = photon_id;
-	}
-} // propagate
-
-} // extern "c"