summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/alpha.h89
-rw-r--r--src/kernel.cu60
-rw-r--r--src/sorting.h103
-rw-r--r--src/transform.cu47
4 files changed, 210 insertions, 89 deletions
diff --git a/src/alpha.h b/src/alpha.h
index e4a3a40..ac75834 100644
--- a/src/alpha.h
+++ b/src/alpha.h
@@ -4,57 +4,14 @@
#include "linalg.h"
#include "intersect.h"
#include "mesh.h"
+#include "sorting.h"
-template <class T>
-__device__ void swap(T &a, T &b)
-{
- T temp = a;
- a = b;
- b = temp;
-}
+#include "stdio.h"
#define ALPHA_DEPTH 10
-struct HitList
-{
- unsigned int size;
- unsigned int indices[ALPHA_DEPTH];
- float distances[ALPHA_DEPTH];
-};
-
-__device__ void add_to_hit_list(const float &distance, const int &index, HitList &h)
-{
- unsigned int i;
- if (h.size >= ALPHA_DEPTH)
- {
- if (distance > h.distances[ALPHA_DEPTH-1])
- return;
-
- i = ALPHA_DEPTH-1;
- }
- else
- {
- i = h.size;
- h.size += 1;
- }
-
- h.indices[i] = index;
- h.distances[i] = distance;
-
- while (i > 0 && h.distances[i-1] > h.distances[i])
- {
- swap(h.distances[i-1], h.distances[i]);
- swap(h.indices[i-1], h.indices[i]);
-
- i -= 1;
- }
-}
-
__device__ int get_color_alpha(const float3 &origin, const float3& direction)
{
- HitList h;
- h.size = 0;
-
float distance;
if (!intersect_node(origin, direction, g_start_node, -1.0f))
@@ -69,6 +26,10 @@ __device__ int get_color_alpha(const float3 &origin, const float3& direction)
unsigned int i;
+ float distances[ALPHA_DEPTH];
+ unsigned int indices[ALPHA_DEPTH];
+ unsigned int n=0;
+
do
{
unsigned int first_child = tex1Dfetch(node_map, *node);
@@ -106,9 +67,26 @@ __device__ int get_color_alpha(const float3 &origin, const float3& direction)
if (intersect_triangle(origin, direction, v0, v1, v2, distance))
{
- add_to_hit_list(distance, i, h);
+ if (n < 1)
+ {
+ distances[0] = distance;
+ indices[0] = i;
+ }
+ else
+ {
+ unsigned long j = searchsorted(n, distances, distance);
+
+ if (j <= ALPHA_DEPTH-1)
+ {
+ insert(ALPHA_DEPTH, distances, j, distance);
+ insert(ALPHA_DEPTH, indices, j, i);
+ }
+ }
+
+ if (n < ALPHA_DEPTH)
+ n++;
}
-
+
} // triangle loop
node--;
@@ -118,16 +96,19 @@ __device__ int get_color_alpha(const float3 &origin, const float3& direction)
} // while loop
while (node != head);
- if (h.size < 1)
+ if (n < 1)
return 0;
float scale = 1.0f;
float fr = 0.0f;
float fg = 0.0f;
float fb = 0.0f;
- for (i=0; i < h.size; i++)
+ unsigned int index;
+ for (i=0; i < n; i++)
{
- uint4 triangle_data = g_triangles[h.indices[i]];
+ index = indices[i];
+
+ uint4 triangle_data = g_triangles[index];
float3 v0 = g_vertices[triangle_data.x];
float3 v1 = g_vertices[triangle_data.y];
@@ -138,11 +119,11 @@ __device__ int get_color_alpha(const float3 &origin, const float3& direction)
if (cos_theta < 0.0f)
cos_theta = dot(-normalize(cross(v1-v0,v2-v1)),-direction);
- unsigned int r0 = 0xff & (g_colors[h.indices[i]] >> 16);
- unsigned int g0 = 0xff & (g_colors[h.indices[i]] >> 8);
- unsigned int b0 = 0xff & g_colors[h.indices[i]];
+ unsigned int r0 = 0xff & (g_colors[index] >> 16);
+ unsigned int g0 = 0xff & (g_colors[index] >> 8);
+ unsigned int b0 = 0xff & g_colors[index];
- float alpha = (255 - (0xff & (g_colors[h.indices[i]] >> 24)))/255.0f;
+ float alpha = (255 - (0xff & (g_colors[index] >> 24)))/255.0f;
fr += r0*scale*cos_theta*alpha;
fg += g0*scale*cos_theta*alpha;
diff --git a/src/kernel.cu b/src/kernel.cu
index 307412e..4418305 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -62,41 +62,6 @@ __device__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max
extern "C"
{
-/* Translate the points `a` by the vector `v` */
-__global__ void translate(int nthreads, float3 *a, float3 v)
-{
- int id = blockIdx.x*blockDim.x + threadIdx.x;
-
- if (id >= nthreads)
- return;
-
- a[id] += v;
-}
-
-/* Rotate the points `a` through an angle `phi` counter-clockwise about the
- axis `axis` (when looking towards +infinity). */
-__global__ void rotate(int nthreads, float3 *a, float phi, float3 axis)
-{
- int id = blockIdx.x*blockDim.x + threadIdx.x;
-
- if (id >= nthreads)
- return;
-
- a[id] = rotate(a[id], phi, axis);
-}
-
-__global__ void rotate_around_point(int nthreads, float3 *a, float phi, float3 axis, float3 point)
-{
- int id = blockIdx.x*blockDim.x + threadIdx.x;
-
- if (id >= nthreads)
- return;
-
- a[id] -= point;
- a[id] = rotate(a[id], phi, axis);
- a[id] += point;
-}
-
__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;
@@ -239,6 +204,28 @@ __global__ void process_image(int nthreads, float3 *image, int *pixels, int nima
} // 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;
@@ -288,6 +275,9 @@ __global__ void ray_trace_alpha(int nthreads, float3 *positions, float3 *directi
float3 direction = directions[id];
direction /= norm(direction);
+ bool hit;
+ float distance;
+
pixels[id] = get_color_alpha(position, direction);
} // ray_trace
diff --git a/src/sorting.h b/src/sorting.h
new file mode 100644
index 0000000..ec16bbe
--- /dev/null
+++ b/src/sorting.h
@@ -0,0 +1,103 @@
+template <class T>
+__device__ void swap(T &a, T &b)
+{
+ T tmp = a;
+ a = b;
+ b = tmp;
+}
+
+template <class T>
+__device__ void reverse(int n, T *a)
+{
+ for (int i=0; i < n/2; i++)
+ swap(a[i],a[n-1-i]);
+}
+
+template <class T>
+__device__ void piksrt(int n, T *arr)
+{
+ int i,j;
+ T a;
+
+ for (j=1; j < n; j++)
+ {
+ a = arr[j];
+ i = j-1;
+ while (i >= 0 && arr[i] > a)
+ {
+ arr[i+1] = arr[i];
+ i--;
+ }
+ arr[i+1] = a;
+ }
+}
+
+template <class T, class U>
+__device__ void piksrt2(int n, T *arr, U *brr)
+{
+ int i,j;
+ T a;
+ U b;
+
+ for (j=1; j < n; j++)
+ {
+ a = arr[j];
+ b = brr[j];
+ i = j-1;
+ while (i >= 0 && arr[i] > a)
+ {
+ arr[i+1] = arr[i];
+ brr[i+1] = brr[i];
+ i--;
+ }
+ arr[i+1] = a;
+ brr[i+1] = b;
+ }
+}
+
+template <class T>
+__device__ unsigned long searchsorted(unsigned long n, T *arr, const T &x)
+{
+ unsigned long ju,jm,jl;
+ int ascnd;
+
+ jl = 0;
+ ju = n;
+
+ ascnd = (arr[n-1] > arr[0]);
+
+ while (ju-jl > 1)
+ {
+ jm = (ju+jl) >> 1;
+
+ if (x >= arr[jm] == ascnd)
+ jl = jm;
+ else
+ ju = jm;
+ }
+
+ if (x <= arr[0])
+ return 0;
+ else if (x == arr[n-1])
+ return n-1;
+ else
+ return ju;
+}
+
+template <class T>
+__device__ void insert(unsigned long n, T *arr, unsigned long i, const T &x)
+{
+ unsigned long j;
+ for (j=n-1; j > i; j--)
+ arr[j] = arr[j-1];
+ arr[i] = x;
+}
+
+template <class T>
+__device__ void add_sorted(unsigned long n, T *arr, const T &x)
+{
+ unsigned long i = searchsorted(n, arr, x);
+
+ if (i < n)
+ insert(n, arr, i, x);
+}
diff --git a/src/transform.cu b/src/transform.cu
new file mode 100644
index 0000000..57bd509
--- /dev/null
+++ b/src/transform.cu
@@ -0,0 +1,47 @@
+//-*-c-*-
+
+#include "linalg.h"
+#include "rotate.h"
+
+extern "C"
+{
+
+/* Translate the points `a` by the vector `v` */
+__global__ void translate(int nthreads, float3 *a, float3 v)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ a[id] += v;
+}
+
+/* Rotate the points `a` through an angle `phi` counter-clockwise about the
+ axis `axis` (when looking towards +infinity). */
+__global__ void rotate(int nthreads, float3 *a, float phi, float3 axis)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ a[id] = rotate(a[id], phi, axis);
+}
+
+/* Rotate the points `a` through an angle `phi` counter-clockwise
+ (when looking towards +infinity along `axis`) about the axis defined
+ by the point `point` and the vector `axis` . */
+__global__ void rotate_around_point(int nthreads, float3 *a, float phi, float3 axis, float3 point)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ a[id] -= point;
+ a[id] = rotate(a[id], phi, axis);
+ a[id] += point;
+}
+
+} // extern "c"