summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/__init__.py6
-rw-r--r--src/daq.cu80
-rw-r--r--src/intersect.h153
-rw-r--r--src/kernel.cu459
-rw-r--r--src/linalg.h121
-rw-r--r--src/materials.h66
-rw-r--r--src/matrix.h223
-rw-r--r--src/physical_constants.h8
-rw-r--r--src/random.h22
-rw-r--r--src/rotate.h26
10 files changed, 1164 insertions, 0 deletions
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000..865d612
--- /dev/null
+++ b/src/__init__.py
@@ -0,0 +1,6 @@
+import os
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+
+kernel = open(dir + '/kernel.cu').read()
+daq = open(dir + '/daq.cu').read()
diff --git a/src/daq.cu b/src/daq.cu
new file mode 100644
index 0000000..c79401c
--- /dev/null
+++ b/src/daq.cu
@@ -0,0 +1,80 @@
+// -*-c++-*-
+#include <curand_kernel.h>
+
+__device__ unsigned int float_to_sortable_int(float f)
+{
+ return __float_as_int(f);
+ //int i = __float_as_int(f);
+ //unsigned int mask = -(int)(i >> 31) | 0x80000000;
+ //return i ^ mask;
+}
+
+__device__ float sortable_int_to_float(unsigned int i)
+{
+ return __int_as_float(i);
+ //unsigned int mask = ((i >> 31) - 1) | 0x80000000;
+ //return __int_as_float(i ^ mask);
+}
+
+
+__device__ curandState daq_rng_states[100000];
+
+extern "C" {
+
+ __global__ void init_daq_rng(int nthreads,
+ unsigned long long seed, unsigned long long offset)
+ {
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ curand_init(seed, id, offset, daq_rng_states+id);
+ }
+
+ __global__ void reset_earliest_time_int(float maxtime,
+ int ntime_ints, unsigned int *time_ints)
+ {
+ int id = threadIdx.x + blockDim.x * blockIdx.x;
+ if (id < ntime_ints) {
+ unsigned int maxtime_int = float_to_sortable_int(maxtime);
+ time_ints[id] = maxtime_int;
+ }
+ }
+
+ __global__ void run_daq(int detection_state, float time_rms,
+ int nphotons, float *photon_times, int *photon_states,
+ int *last_hit_triangles, int *solid_map,
+ int nsolids, unsigned int *earliest_time_int)
+ {
+ int id = threadIdx.x + blockDim.x * blockIdx.x;
+
+ curandState_t rng = daq_rng_states[id];
+
+ if (id < nphotons) {
+ int triangle_id = last_hit_triangles[id];
+
+ if (triangle_id > -1) {
+ int solid_id = solid_map[triangle_id];
+ int state = photon_states[id];
+ float time = photon_times[id];// + curand_normal(&rng) * time_rms;
+ unsigned int time_int = float_to_sortable_int(time);
+ if (solid_id < nsolids && state == detection_state)
+ atomicMin(earliest_time_int + solid_id, time_int);
+ }
+ }
+
+ daq_rng_states[id] = rng;
+ }
+
+ __global__ void convert_sortable_int_to_float(int n,
+ unsigned int *sortable_ints,
+ float *float_output)
+ {
+ int id = threadIdx.x + blockDim.x * blockIdx.x;
+
+ if (id < n)
+ float_output[id] = sortable_int_to_float(sortable_ints[id]);
+ }
+
+} // extern "C"
diff --git a/src/intersect.h b/src/intersect.h
new file mode 100644
index 0000000..0b713c8
--- /dev/null
+++ b/src/intersect.h
@@ -0,0 +1,153 @@
+//-*-c-*-
+
+#ifndef __INTERSECT_H__
+#define __INTERSECT_H__
+
+#include <math_constants.h>
+
+#include "linalg.h"
+#include "matrix.h"
+#include "rotate.h"
+
+/* Test the intersection between a ray starting from `origin` traveling in the
+ direction `direction` and a triangle defined by the vertices `v0`, `v1`, and
+ `v2`. If the ray intersects the triangle, set `distance` to the distance
+ between `origin` and the intersection and return true, else return false.
+
+ `direction` must be normalized. */
+__device__ bool intersect_triangle(const float3 &origin, const float3 &direction, const float3 &v0, const float3 &v1, const float3 &v2, float &distance)
+{
+ Matrix m = make_matrix(v1-v0, v2-v0, -direction);
+
+ float determinant = det(m);
+
+ if (determinant == 0.0f)
+ return false;
+
+ float3 b = origin-v0;
+
+ float u1 = ((m.a11*m.a22 - m.a12*m.a21)*b.x +
+ (m.a02*m.a21 - m.a01*m.a22)*b.y +
+ (m.a01*m.a12 - m.a02*m.a11)*b.z)/determinant;
+
+ if (u1 < 0.0f || u1 > 1.0f)
+ return false;
+
+ float u2 = ((m.a12*m.a20 - m.a10*m.a22)*b.x +
+ (m.a00*m.a22 - m.a02*m.a20)*b.y +
+ (m.a02*m.a10 - m.a00*m.a12)*b.z)/determinant;
+
+ if (u2 < 0.0f || u2 > 1.0f)
+ return false;
+
+ float u3 = ((m.a10*m.a21 - m.a11*m.a20)*b.x +
+ (m.a01*m.a20 - m.a00*m.a21)*b.y +
+ (m.a00*m.a11 - m.a01*m.a10)*b.z)/determinant;
+
+ if (u3 <= 0.0f || (1.0f-u1-u2) < 0.0f)
+ return false;
+
+ distance = u3;
+
+ return true;
+}
+
+/* Return the 32 bit color associated with the intersection between a ray
+ starting from `origin` traveling in the direction `direction` and the
+ plane defined by the points `v0`, `v1`, and `v2` using the cosine of the
+ angle between the ray and the plane normal to determine the brightness.
+
+ `direction` must be normalized. */
+__device__ int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2, const int base_color=0xFFFFFFFF)
+{
+ float scale = dot(normalize(cross(v1-v0,v2-v1)),-direction);
+
+ unsigned int r = 0xFF & (base_color >> 16);
+ unsigned int g = 0xFF & (base_color >> 8);
+ unsigned int b = 0xFF & base_color;
+
+ if (scale < 0.0f)
+ scale = dot(-normalize(cross(v1-v0,v2-v0)),-direction);
+
+ r = floorf(r*scale);
+ g = floorf(g*scale);
+ b = floorf(b*scale);
+
+ return r << 16 | g << 8 | b;
+}
+
+/* Test the intersection between a ray starting from `origin` traveling in the
+ direction `direction` and the axis-aligned box defined by the opposite
+ vertices `lower_bound` and `upper_bound`. If the ray intersects the box
+ return True, else return False. */
+__device__ bool intersect_box(const float3 &origin, const float3 &direction, const float3 &lower_bound, const float3 &upper_bound)
+{
+ float kmin, kmax, kymin, kymax, kzmin, kzmax;
+
+ float divx = 1.0f/direction.x;
+ if (divx >= 0.0f)
+ {
+ kmin = (lower_bound.x - origin.x)*divx;
+ kmax = (upper_bound.x - origin.x)*divx;
+ }
+ else
+ {
+ kmin = (upper_bound.x - origin.x)*divx;
+ kmax = (lower_bound.x - origin.x)*divx;
+ }
+
+ if (kmax < 0.0f)
+ return false;
+
+ float divy = 1.0f/direction.y;
+ if (divy >= 0.0f)
+ {
+ kymin = (lower_bound.y - origin.y)*divy;
+ kymax = (upper_bound.y - origin.y)*divy;
+ }
+ else
+ {
+ kymin = (upper_bound.y - origin.y)*divy;
+ kymax = (lower_bound.y - origin.y)*divy;
+ }
+
+ if (kymax < 0.0f)
+ return false;
+
+ if (kymin > kmin)
+ kmin = kymin;
+
+ if (kymax < kmax)
+ kmax = kymax;
+
+ if (kmin > kmax)
+ return false;
+
+ float divz = 1.0f/direction.z;
+ if (divz >= 0.0f)
+ {
+ kzmin = (lower_bound.z - origin.z)*divz;
+ kzmax = (upper_bound.z - origin.z)*divz;
+ }
+ else
+ {
+ kzmin = (upper_bound.z - origin.z)*divz;
+ kzmax = (lower_bound.z - origin.z)*divz;
+ }
+
+ if (kzmax < 0.0f)
+ return false;
+
+ if (kzmin > kmin)
+ kmin = kzmin;
+
+ if (kzmax < kmax)
+ kmax = kzmax;
+
+ if (kmin > kmax)
+ return false;
+
+ return true;
+}
+
+#endif
diff --git a/src/kernel.cu b/src/kernel.cu
new file mode 100644
index 0000000..2d579be
--- /dev/null
+++ b/src/kernel.cu
@@ -0,0 +1,459 @@
+//-*-c-*-
+#include <math_constants.h>
+#include <curand_kernel.h>
+
+#include "linalg.h"
+#include "matrix.h"
+#include "rotate.h"
+#include "intersect.h"
+#include "materials.h"
+#include "physical_constants.h"
+#include "random.h"
+
+#define STACK_SIZE 500
+
+enum
+{
+ DEBUG = -2,
+ INIT = -1,
+ NO_HIT,
+ BULK_ABSORB,
+ SURFACE_ABSORB
+};
+
+/* flattened triangle mesh */
+__device__ float4 *vertices;
+__device__ uint4 *triangles;
+
+/* lower/upper bounds for the bounding box associated with each node/leaf */
+texture<float4, 1, cudaReadModeElementType> upper_bounds;
+texture<float4, 1, cudaReadModeElementType> lower_bounds;
+
+/* map to child nodes/triangles and the number of child nodes/triangles */
+texture<unsigned int, 1, cudaReadModeElementType> node_map;
+texture<unsigned int, 1, cudaReadModeElementType> node_length;
+
+__device__ float3 make_float3(const float4 &a)
+{
+ return make_float3(a.x, a.y, a.z);
+}
+
+__device__ int convert(int c)
+{
+ if (c & 0x80)
+ return (0xFFFFFF00 | c);
+ else
+ return c;
+}
+
+/* Test the intersection between a ray starting at `origin` traveling in the
+ direction `direction` and the bounding box around node `i`. If the ray
+ intersects the bounding box return true, else return false. */
+__device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i)
+{
+ float3 lower_bound = make_float3(tex1Dfetch(lower_bounds, i));
+ float3 upper_bound = make_float3(tex1Dfetch(upper_bounds, i));
+
+ return intersect_box(origin, direction, lower_bound, upper_bound);
+}
+
+/* Find the intersection between a ray starting at `origin` traveling in the
+ 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, float &min_distance, int last_hit_triangle = -1)
+{
+ int triangle_index = -1;
+
+ float distance;
+
+ if (!intersect_node(origin, direction, start_node))
+ return -1;
+
+ int stack[STACK_SIZE];
+
+ int *head = &stack[0];
+ int *node = &stack[1];
+ int *tail = &stack[STACK_SIZE-1];
+ *node = start_node;
+
+ int i;
+
+ do
+ {
+ int index = tex1Dfetch(node_map, *node);
+ int length = tex1Dfetch(node_length, *node);
+
+ if (*node >= first_node)
+ {
+ for (i=0; i < length; i++)
+ {
+ if (intersect_node(origin, direction, index+i))
+ {
+ //*node++ = index+i;
+ *node = index+i;
+ node++;
+ }
+ }
+
+ if (node == head)
+ break;
+
+ node--;
+ }
+ else // node is a leaf
+ {
+ for (i=0; i < length; i++)
+ {
+ if (last_hit_triangle == index+i)
+ continue;
+
+ uint4 triangle_data = triangles[index+i];
+
+ float3 v0 = make_float3(vertices[triangle_data.x]);
+ float3 v1 = make_float3(vertices[triangle_data.y]);
+ float3 v2 = make_float3(vertices[triangle_data.z]);
+
+ if (intersect_triangle(origin, direction, v0, v1, v2, distance))
+ {
+ if (triangle_index == -1)
+ {
+ triangle_index = index + i;
+ min_distance = distance;
+ continue;
+ }
+
+ if (distance < min_distance)
+ {
+ triangle_index = index + i;
+ min_distance = distance;
+ }
+ }
+ } // triangle loop
+
+ node--;
+
+ } // node is a leaf
+
+ } // while loop
+ while (node != head);
+
+ return triangle_index;
+}
+
+__device__ curandState rng_states[100000];
+
+extern "C"
+{
+
+ __global__ void set_pointer(uint4 *triangle_ptr, float4 *vertex_ptr)
+{
+ triangles = triangle_ptr;
+ vertices = vertex_ptr;
+}
+
+/* Initialize random number states */
+__global__ void init_rng(int nthreads, unsigned long long seed, unsigned long long offset)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ curand_init(seed, id, offset, rng_states+id);
+}
+
+/* Translate `points` by the vector `v` */
+__global__ void translate(int nthreads, float3 *points, float3 v)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ points[id] += v;
+}
+
+/* Rotate `points` through an angle `phi` counter-clockwise about the
+ axis `axis` (when looking towards +infinity). */
+__global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ points[id] = rotate(points[id], phi, axis);
+}
+
+/* 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(int nthreads, float3 *positions, float3 *directions, int start_node, int first_node, 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, start_node, first_node, distance);
+
+ if (triangle_index == -1)
+ {
+ pixels[id] = 0x000000;
+ }
+ else
+ {
+ uint4 triangle_data = triangles[triangle_index];
+
+ float3 v0 = make_float3(vertices[triangle_data.x]);
+ float3 v1 = make_float3(vertices[triangle_data.y]);
+ float3 v2 = make_float3(vertices[triangle_data.z]);
+
+ pixels[id] = get_color(direction, v0, v1, v2, triangle_data.w);
+ }
+
+} // ray_trace
+
+__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 max_steps)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ int state = states[id];
+
+ if (state != INIT)
+ return;
+
+ curandState rng = rng_states[id];
+ float3 position = positions[id];
+ float3 direction = directions[id];
+ float3 polarization = polarizations[id];
+ float wavelength = wavelengths[id];
+ float time = times[id];
+ int last_hit_triangle = last_hit_triangles[id];
+
+ int steps = 0;
+ while (steps < max_steps)
+ {
+ steps++;
+
+ direction /= norm(direction);
+ polarization /= norm(polarization);
+
+ float distance;
+
+ last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
+
+ if (last_hit_triangle == -1)
+ {
+ state = NO_HIT;
+ break;
+ }
+
+ uint4 triangle_data = triangles[last_hit_triangle];
+
+ float3 v0 = make_float3(vertices[triangle_data.x]);
+ float3 v1 = make_float3(vertices[triangle_data.y]);
+ float3 v2 = make_float3(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 = -absorption_length*logf(curand_uniform(&rng));
+
+ if (absorption_distance <= scattering_distance)
+ {
+ if (absorption_distance <= distance)
+ {
+ time += absorption_distance/(SPEED_OF_LIGHT/refractive_index1);
+ position += absorption_distance*direction;
+ state = BULK_ABSORB;
+
+ last_hit_triangle = -1;
+
+ break;
+ } // 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 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 uniform_sample = curand_uniform(&rng);
+
+ if (uniform_sample < absorption)
+ {
+ // absorb
+ state = SURFACE_ABSORB;
+ break;
+ }
+ else if (uniform_sample < absorption + reflection_diffuse)
+ {
+ // diffusely reflect
+ direction = uniform_sphere(&rng);
+
+ if (dot(direction, surface_normal) < 0.0f)
+ direction = -direction;
+
+ // randomize polarization ?
+ polarization = cross(uniform_sphere(&rng), direction);
+ polarization /= norm(polarization);
+
+ continue;
+ }
+ else
+ {
+ // specularly reflect
+
+ direction = rotate(surface_normal, incident_angle, p);
+
+ 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);
+
+ continue;
+ } // photon polarization parallel to surface
+
+ } // while(nsteps < max_steps)
+
+ states[id] = state;
+ rng_states[id] = rng;
+ positions[id] = position;
+ directions[id] = direction;
+ polarizations[id] = polarization;
+ wavelengths[id] = wavelength;
+ times[id] = time;
+ last_hit_triangles[id] = last_hit_triangle;
+
+} // propagate
+
+} // extern "c"
diff --git a/src/linalg.h b/src/linalg.h
new file mode 100644
index 0000000..fe6b1b2
--- /dev/null
+++ b/src/linalg.h
@@ -0,0 +1,121 @@
+#ifndef __LINALG_H__
+#define __LINALG_H__
+
+__device__ __host__ float3 operator- (const float3 &a)
+{
+ return make_float3(-a.x, -a.y, -a.z);
+}
+
+__device__ __host__ float3 operator+ (const float3 &a, const float3 &b)
+{
+ return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
+}
+
+__device__ __host__ void operator+= (float3 &a, const float3 &b)
+{
+ a.x += b.x;
+ a.y += b.y;
+ a.z += b.z;
+}
+
+__device__ __host__ float3 operator- (const float3 &a, const float3 &b)
+{
+ return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
+}
+
+__device__ __host__ void operator-= (float3 &a, const float3 &b)
+{
+ a.x -= b.x;
+ a.y -= b.y;
+ a.z -= b.z;
+}
+
+__device__ __host__ float3 operator+ (const float3 &a, const float &c)
+{
+ return make_float3(a.x+c, a.y+c, a.z+c);
+}
+
+__device__ __host__ void operator+= (float3 &a, const float &c)
+{
+ a.x += c;
+ a.y += c;
+ a.z += c;
+}
+
+__device__ __host__ float3 operator+ (const float &c, const float3 &a)
+{
+ return make_float3(c+a.x, c+a.y, c+a.z);
+}
+
+__device__ __host__ float3 operator- (const float3 &a, const float &c)
+{
+ return make_float3(a.x-c, a.y-c, a.z-c);
+}
+
+__device__ __host__ void operator-= (float3 &a, const float &c)
+{
+ a.x -= c;
+ a.y -= c;
+ a.z -= c;
+}
+
+__device__ __host__ float3 operator- (const float &c, const float3& a)
+{
+ return make_float3(c-a.x, c-a.y, c-a.z);
+}
+
+__device__ __host__ float3 operator* (const float3 &a, const float &c)
+{
+ return make_float3(a.x*c, a.y*c, a.z*c);
+}
+
+__device__ __host__ void operator*= (float3 &a, const float &c)
+{
+ a.x *= c;
+ a.y *= c;
+ a.z *= c;
+}
+
+__device__ __host__ float3 operator* (const float &c, const float3& a)
+{
+ return make_float3(c*a.x, c*a.y, c*a.z);
+}
+
+__device__ __host__ float3 operator/ (const float3 &a, const float &c)
+{
+ return make_float3(a.x/c, a.y/c, a.z/c);
+}
+
+__device__ __host__ void operator/= (float3 &a, const float &c)
+{
+ a.x /= c;
+ a.y /= c;
+ a.z /= c;
+}
+
+__device__ __host__ float3 operator/ (const float &c, const float3 &a)
+{
+ return make_float3(c/a.x, c/a.y, c/a.z);
+}
+
+__device__ __host__ float dot(const float3 &a, const float3 &b)
+{
+ return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+__device__ __host__ float3 cross(const float3 &a, const float3 &b)
+{
+ return make_float3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
+}
+
+__device__ __host__ float norm(const float3 &a)
+{
+ return sqrtf(dot(a,a));
+}
+
+__device__ __host__ float3 normalize(const float3 &a)
+{
+ return a/norm(a);
+}
+
+#endif
diff --git a/src/materials.h b/src/materials.h
new file mode 100644
index 0000000..77c9f43
--- /dev/null
+++ b/src/materials.h
@@ -0,0 +1,66 @@
+#ifndef __MATERIALS_H__
+#define __MATERIALS_H__
+
+__device__ float min_wavelength;
+__device__ float max_wavelength;
+__device__ float wavelength_step;
+__device__ unsigned int wavelength_size;
+
+struct Material
+{
+ float *refractive_index;
+ float *absorption_length;
+ float *scattering_length;
+};
+
+struct Surface
+{
+ float *absorption;
+ float *reflection_diffuse;
+ float *reflection_specular;
+};
+
+__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];
+
+ int jl = (x-min_wavelength)/wavelength_step;
+
+ return fp[jl] + (x-(min_wavelength + jl*wavelength_step))*(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;
+ materials[material_index].absorption_length = absorption_length;
+ materials[material_index].scattering_length = scattering_length;
+}
+
+__global__ void set_surface(int surface_index, float *absorption, float *reflection_diffuse, float *reflection_specular)
+{
+ surfaces[surface_index].absorption = absorption;
+ surfaces[surface_index].reflection_diffuse = reflection_diffuse;
+ surfaces[surface_index].reflection_specular = reflection_specular;
+}
+
+} // extern "c"
+
+#endif
diff --git a/src/matrix.h b/src/matrix.h
new file mode 100644
index 0000000..a3e480b
--- /dev/null
+++ b/src/matrix.h
@@ -0,0 +1,223 @@
+#ifndef __MATRIX_H__
+#define __MATRIX_H__
+
+struct Matrix
+{
+ float a00, a01, a02, a10, a11, a12, a20, a21, a22;
+};
+
+__device__ __host__ Matrix make_matrix(float a00, float a01, float a02,
+ float a10, float a11, float a12,
+ float a20, float a21, float a22)
+{
+ Matrix m = {a00, a01, a02, a10, a11, a12, a20, a21, a22};
+ return m;
+}
+
+__device__ __host__ Matrix make_matrix(const float3 &u1, const float3 &u2, const float3 &u3)
+{
+ Matrix m = {u1.x, u2.x, u3.x, u1.y, u2.y, u3.y, u1.z, u2.z, u3.z};
+ return m;
+}
+
+__device__ __host__ Matrix operator- (const Matrix &m)
+{
+ return make_matrix(-m.a00, -m.a01, -m.a02,
+ -m.a10, -m.a11, -m.a12,
+ -m.a20, -m.a21, -m.a22);
+}
+
+__device__ __host__ float3 operator* (const Matrix &m, const float3 &a)
+{
+ return make_float3(m.a00*a.x + m.a01*a.y + m.a02*a.z,
+ m.a10*a.x + m.a11*a.y + m.a12*a.z,
+ m.a20*a.x + m.a21*a.y + m.a22*a.z);
+}
+
+__device__ __host__ Matrix operator+ (const Matrix &m, const Matrix &n)
+{
+ return make_matrix(m.a00+n.a00, m.a01+n.a01, m.a02+n.a02,
+ m.a10+n.a10, m.a11+n.a11, m.a12+n.a12,
+ m.a20+n.a20, m.a21+n.a21, m.a22+n.a22);
+}
+
+__device__ __host__ void operator+= (Matrix &m, const Matrix &n)
+{
+ m.a00 += n.a00;
+ m.a01 += n.a01;
+ m.a02 += n.a02;
+ m.a10 += n.a10;
+ m.a11 += n.a11;
+ m.a12 += n.a12;
+ m.a20 += n.a20;
+ m.a21 += n.a21;
+ m.a22 += n.a22;
+}
+
+__device__ __host__ Matrix operator- (const Matrix &m, const Matrix &n)
+{
+ return make_matrix(m.a00-n.a00, m.a01-n.a01, m.a02-n.a02,
+ m.a10-n.a10, m.a11-n.a11, m.a12-n.a12,
+ m.a20-n.a20, m.a21-n.a21, m.a22-n.a22);
+}
+
+__device__ __host__ void operator-= (Matrix &m, const Matrix& n)
+{
+ m.a00 -= n.a00;
+ m.a01 -= n.a01;
+ m.a02 -= n.a02;
+ m.a10 -= n.a10;
+ m.a11 -= n.a11;
+ m.a12 -= n.a12;
+ m.a20 -= n.a20;
+ m.a21 -= n.a21;
+ m.a22 -= n.a22;
+}
+
+__device__ __host__ Matrix operator* (const Matrix &m, const Matrix &n)
+{
+ return make_matrix(m.a00*n.a00 + m.a01*n.a10 + m.a02*n.a20,
+ m.a00*n.a01 + m.a01*n.a11 + m.a02*n.a21,
+ m.a00*n.a02 + m.a01*n.a12 + m.a02*n.a22,
+ m.a10*n.a00 + m.a11*n.a10 + m.a12*n.a20,
+ m.a10*n.a01 + m.a11*n.a11 + m.a12*n.a21,
+ m.a10*n.a02 + m.a11*n.a12 + m.a12*n.a22,
+ m.a20*n.a00 + m.a21*n.a10 + m.a22*n.a20,
+ m.a20*n.a01 + m.a21*n.a11 + m.a22*n.a21,
+ m.a20*n.a02 + m.a21*n.a12 + m.a22*n.a22);
+}
+
+__device__ __host__ Matrix operator+ (const Matrix &m, const float &c)
+{
+ return make_matrix(m.a00+c, m.a01+c, m.a02+c,
+ m.a10+c, m.a11+c, m.a12+c,
+ m.a20+c, m.a21+c, m.a22+c);
+}
+
+__device__ __host__ void operator+= (Matrix &m, const float &c)
+{
+ m.a00 += c;
+ m.a01 += c;
+ m.a02 += c;
+ m.a10 += c;
+ m.a11 += c;
+ m.a12 += c;
+ m.a20 += c;
+ m.a21 += c;
+ m.a22 += c;
+}
+
+__device__ __host__ Matrix operator+ (const float &c, const Matrix &m)
+{
+ return make_matrix(c+m.a00, c+m.a01, c+m.a02,
+ c+m.a10, c+m.a11, c+m.a12,
+ c+m.a20, c+m.a21, c+m.a22);
+}
+
+__device__ __host__ Matrix operator- (const Matrix &m, const float &c)
+{
+ return make_matrix(m.a00-c, m.a01-c, m.a02-c,
+ m.a10-c, m.a11-c, m.a12-c,
+ m.a20-c, m.a21-c, m.a22-c);
+}
+
+__device__ __host__ void operator-= (Matrix &m, const float &c)
+{
+ m.a00 -= c;
+ m.a01 -= c;
+ m.a02 -= c;
+ m.a10 -= c;
+ m.a11 -= c;
+ m.a12 -= c;
+ m.a20 -= c;
+ m.a21 -= c;
+ m.a22 -= c;
+}
+
+__device__ __host__ Matrix operator- (const float &c, const Matrix &m)
+{
+ return make_matrix(c-m.a00, c-m.a01, c-m.a02,
+ c-m.a10, c-m.a11, c-m.a12,
+ c-m.a20, c-m.a21, c-m.a22);
+}
+
+__device__ __host__ Matrix operator* (const Matrix &m, const float &c)
+{
+ return make_matrix(m.a00*c, m.a01*c, m.a02*c,
+ m.a10*c, m.a11*c, m.a12*c,
+ m.a20*c, m.a21*c, m.a22*c);
+}
+
+__device__ __host__ void operator*= (Matrix &m, const float &c)
+{
+ m.a00 *= c;
+ m.a01 *= c;
+ m.a02 *= c;
+ m.a10 *= c;
+ m.a11 *= c;
+ m.a12 *= c;
+ m.a20 *= c;
+ m.a21 *= c;
+ m.a22 *= c;
+}
+
+__device__ __host__ Matrix operator* (const float &c, const Matrix &m)
+{
+ return make_matrix(c*m.a00, c*m.a01, c*m.a02,
+ c*m.a10, c*m.a11, c*m.a12,
+ c*m.a20, c*m.a21, c*m.a22);
+}
+
+__device__ __host__ Matrix operator/ (const Matrix &m, const float &c)
+{
+ return make_matrix(m.a00/c, m.a01/c, m.a02/c,
+ m.a10/c, m.a11/c, m.a12/c,
+ m.a20/c, m.a21/c, m.a22/c);
+}
+
+__device__ __host__ void operator/= (Matrix &m, const float &c)
+{
+ m.a00 /= c;
+ m.a01 /= c;
+ m.a02 /= c;
+ m.a10 /= c;
+ m.a11 /= c;
+ m.a12 /= c;
+ m.a20 /= c;
+ m.a21 /= c;
+ m.a22 /= c;
+}
+
+__device__ __host__ Matrix operator/ (const float &c, const Matrix &m)
+{
+ return make_matrix(c/m.a00, c/m.a01, c/m.a02,
+ c/m.a10, c/m.a11, c/m.a12,
+ c/m.a20, c/m.a21, c/m.a22);
+}
+
+__device__ __host__ float det(const Matrix &m)
+{
+ return m.a00*(m.a11*m.a22 - m.a12*m.a21) - m.a10*(m.a01*m.a22 - m.a02*m.a21) + m.a20*(m.a01*m.a12 - m.a02*m.a11);
+}
+
+__device__ __host__ Matrix inv(const Matrix &m)
+{
+ return make_matrix(m.a11*m.a22 - m.a12*m.a21,
+ m.a02*m.a21 - m.a01*m.a22,
+ m.a01*m.a12 - m.a02*m.a11,
+ m.a12*m.a20 - m.a10*m.a22,
+ m.a00*m.a22 - m.a02*m.a20,
+ m.a02*m.a10 - m.a00*m.a12,
+ m.a10*m.a21 - m.a11*m.a20,
+ m.a01*m.a20 - m.a00*m.a21,
+ m.a00*m.a11 - m.a01*m.a10)/det(m);
+}
+
+__device__ __host__ Matrix outer(const float3 &a, const float3 &b)
+{
+ return make_matrix(a.x*b.x, a.x*b.y, a.x*b.z,
+ a.y*b.x, a.y*b.y, a.y*b.z,
+ a.z*b.x, a.z*b.y, a.z*b.z);
+}
+
+#endif
diff --git a/src/physical_constants.h b/src/physical_constants.h
new file mode 100644
index 0000000..d591e93
--- /dev/null
+++ b/src/physical_constants.h
@@ -0,0 +1,8 @@
+#ifndef __PHYSICAL_CONSTANTS_H__
+#define __PHYSICAL_CONSTANTS_H__
+
+#define SPEED_OF_LIGHT 2.99792458e8
+#define PI 3.141592653589793
+#define EPSILON 1.0e-3
+
+#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
diff --git a/src/rotate.h b/src/rotate.h
new file mode 100644
index 0000000..7fc6f7e
--- /dev/null
+++ b/src/rotate.h
@@ -0,0 +1,26 @@
+#ifndef __ROTATE_H__
+#define __ROTATE_H__
+
+#include "linalg.h"
+#include "matrix.h"
+
+__device__ const Matrix IDENTITY_MATRIX = {1,0,0,0,1,0,0,0,1};
+
+__device__ __host__ Matrix make_rotation_matrix(float phi, const float3 &n)
+{
+ /* rotate points counterclockwise, when looking towards +infinity,
+ through an angle `phi` about the axis `n`. */
+
+ float cos_phi = cosf(phi);
+ float sin_phi = sinf(phi);
+
+ return IDENTITY_MATRIX*cos_phi + (1-cos_phi)*outer(n,n) +
+ sin_phi*make_matrix(0,n.z,-n.y,-n.z,0,n.x,n.y,-n.x,0);
+}
+
+__device__ __host__ float3 rotate(const float3 &a, float phi, const float3 &n)
+{
+ return make_rotation_matrix(phi,n)*a;
+}
+
+#endif
generated by cgit (git 2.34.1) at 2025-05-10 04:00:27 +0000