From 9306f888fea903accf827870a122a2f6f76e472e Mon Sep 17 00:00:00 2001 From: Anthony LaTorre Date: Wed, 18 May 2011 11:29:26 -0400 Subject: added some more documentation and a more accurate miniature version of lbne --- src/intersect.cu | 300 ------------------------------------------------------- src/intersect.h | 141 ++++++++++++++++++++++++++ src/kernel.cu | 193 +++++++++++++++++++++++++++++++++++ 3 files changed, 334 insertions(+), 300 deletions(-) delete mode 100644 src/intersect.cu create mode 100644 src/intersect.h create mode 100644 src/kernel.cu (limited to 'src') diff --git a/src/intersect.cu b/src/intersect.cu deleted file mode 100644 index 1402aa1..0000000 --- a/src/intersect.cu +++ /dev/null @@ -1,300 +0,0 @@ -//-*-c-*- -#include - -#include "linalg.h" -#include "matrix.h" -#include "rotate.h" - -/* flattened triangle mesh */ -texture mesh; - -/* lower/upper bounds for the bounding box associated with each node/leaf */ -texture upper_bound_arr; -texture lower_bound_arr; - -/* map to child nodes/triangles and the number of child nodes/triangles */ -texture child_map_arr; -texture child_len_arr; - -__device__ float3 make_float3(const float4 &a) -{ - return make_float3(a.x, a.y, a.z); -} - -/* 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) - 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) - 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-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) -{ - float scale = 255*dot(normalize(cross(v1-v0,v2-v0)),-direction); - - if (scale < 0.0f) - scale = 255*dot(-normalize(cross(v1-v0,v2-v0)),-direction); - - int rgb = floorf(scale); - - return rgb << 16 | rgb << 8 | rgb; -} - -/* 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; - - if (direction.x >= 0.0f) - { - kmin = (lower_bound.x - origin.x)/direction.x; - kmax = (upper_bound.x - origin.x)/direction.x; - } - else - { - kmin = (upper_bound.x - origin.x)/direction.x; - kmax = (lower_bound.x - origin.x)/direction.x; - } - - if (kmax < kmin) - return false; - - if (direction.y >= 0.0f) - { - kymin = (lower_bound.y - origin.y)/direction.y; - kymax = (upper_bound.y - origin.y)/direction.y; - } - else - { - kymin = (upper_bound.y - origin.y)/direction.y; - kymax = (lower_bound.y - origin.y)/direction.y; - } - - if (kymax < kymin) - return false; - - if (kymin > kmin) - kmin = kymin; - - if (kymax < kmax) - kmax = kymax; - - if (kmin > kmax) - return false; - - if (direction.z >= 0.0f) - { - kzmin = (lower_bound.z - origin.z)/direction.z; - kzmax = (upper_bound.z - origin.z)/direction.z; - } - else - { - kzmin = (upper_bound.z - origin.z)/direction.z; - kzmax = (lower_bound.z - origin.z)/direction.z; - } - - if (kzmax < kzmin) - return false; - - if (kzmin > kmin) - kmin = kzmin; - - if (kzmax < kmax) - kmax = kzmax; - - if (kmin > kmax) - return false; - - if (kmax < 0.0f) - return false; - - return true; -} - -/* 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_bound_arr, i)); - float3 upper_bound = make_float3(tex1Dfetch(upper_bound_arr, i)); - - return intersect_box(origin, direction, lower_bound, upper_bound); -} - -extern "C" -{ - -__global__ void translate(int max_idx, float3 *pt, float3 v) -{ - int idx = blockIdx.x*blockDim.x + threadIdx.x; - - if (idx > max_idx) - return; - - pt[idx] += v; -} - -__global__ void rotate(int max_idx, float3 *pt, float phi, float3 axis) -{ - int idx = blockIdx.x*blockDim.x + threadIdx.x; - - if (idx > max_idx) - return; - - pt[idx] = rotate(pt[idx], phi, axis); -} - -__global__ void intersect_mesh(int max_idx, float3 *origin_arr, float3 *direction_arr, int first_leaf, int *state_arr, int *pixel_arr) -{ - int idx = blockIdx.x*blockDim.x + threadIdx.x; - - if (idx > max_idx) - return; - - float3 origin = origin_arr[idx]; - float3 direction = direction_arr[idx]; - direction /= norm(direction); - - int *pixel = pixel_arr+idx; - int *state = state_arr+idx; - - bool hit = false; - - float distance; - float min_distance; - - if (!intersect_node(origin, direction, 0)) - { - *pixel = 0; - return; - } - - int stack[100]; - - int *head = &stack[0]; - int *node = &stack[1]; - *node = 0; - - int i; - - bool show_leafs = false; - - do - { - int child_map = tex1Dfetch(child_map_arr, *node); - int child_len = tex1Dfetch(child_len_arr, *node); - - if (*node < first_leaf) - { - for (i=0; i < child_len; i++) - if (intersect_node(origin, direction, child_map+i)) - *node++ = child_map+i; - - if (node == head) - break; - - node--; - - } - else if (show_leafs) - { - hit = true; - *pixel = 255; - node--; - } - else // node is a leaf - { - for (i=0; i < child_len; i++) - { - int mesh_idx = 3*(child_map + i); - - float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx)); - float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1)); - float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2)); - - if (intersect_triangle(origin, direction, v0, v1, v2, distance)) - { - if (!hit) - { - *pixel = get_color(direction, v0, v1, v2); - *state = child_map + i; - - min_distance = distance; - - hit = true; - continue; - } - - if (distance < min_distance) - { - *pixel = get_color(direction, v0, v1, v2); - *state = child_map + i; - - min_distance = distance; - } - } - - } // triangle loop - - node--; - - } // node is a leaf - - } // while loop - while (node != head); - - if (!hit) - { - *state = -1; - *pixel = 0; - } - -} // intersect mesh - -} // extern "c" diff --git a/src/intersect.h b/src/intersect.h new file mode 100644 index 0000000..b984612 --- /dev/null +++ b/src/intersect.h @@ -0,0 +1,141 @@ +//-*-c-*- +#include + +#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) + 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) + 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-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) +{ + float scale = 255*dot(normalize(cross(v1-v0,v2-v0)),-direction); + + if (scale < 0.0f) + scale = 255*dot(-normalize(cross(v1-v0,v2-v0)),-direction); + + int rgb = floorf(scale); + + return rgb << 16 | rgb << 8 | rgb; +} + +/* 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; + + if (direction.x >= 0.0f) + { + kmin = (lower_bound.x - origin.x)/direction.x; + kmax = (upper_bound.x - origin.x)/direction.x; + } + else + { + kmin = (upper_bound.x - origin.x)/direction.x; + kmax = (lower_bound.x - origin.x)/direction.x; + } + + if (kmax < kmin) + return false; + + if (direction.y >= 0.0f) + { + kymin = (lower_bound.y - origin.y)/direction.y; + kymax = (upper_bound.y - origin.y)/direction.y; + } + else + { + kymin = (upper_bound.y - origin.y)/direction.y; + kymax = (lower_bound.y - origin.y)/direction.y; + } + + if (kymax < kymin) + return false; + + if (kymin > kmin) + kmin = kymin; + + if (kymax < kmax) + kmax = kymax; + + if (kmin > kmax) + return false; + + if (direction.z >= 0.0f) + { + kzmin = (lower_bound.z - origin.z)/direction.z; + kzmax = (upper_bound.z - origin.z)/direction.z; + } + else + { + kzmin = (upper_bound.z - origin.z)/direction.z; + kzmax = (lower_bound.z - origin.z)/direction.z; + } + + if (kzmax < kzmin) + return false; + + if (kzmin > kmin) + kmin = kzmin; + + if (kzmax < kmax) + kmax = kzmax; + + if (kmin > kmax) + return false; + + if (kmax < 0.0f) + return false; + + return true; +} diff --git a/src/kernel.cu b/src/kernel.cu new file mode 100644 index 0000000..ed53032 --- /dev/null +++ b/src/kernel.cu @@ -0,0 +1,193 @@ +//-*-c-*- +#include + +#include "linalg.h" +#include "matrix.h" +#include "rotate.h" +#include "intersect.h" + +#define STACK_SIZE 100 + +/* flattened triangle mesh */ +texture mesh; + +/* lower/upper bounds for the bounding box associated with each node/leaf */ +texture upper_bounds; +texture lower_bounds; + +/* map to child nodes/triangles and the number of child nodes/triangles */ +texture child_map_arr; +texture child_len_arr; + +__device__ float3 make_float3(const float4 &a) +{ + return make_float3(a.x, a.y, a.z); +} + +/* 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 first_leaf) +{ + int triangle_idx = -1; + + float distance; + float min_distance; + + if (!intersect_node(origin, direction, 0)) + return -1; + + int stack[STACK_SIZE]; + + int *head = &stack[0]; + int *node = &stack[1]; + int *tail = &stack[STACK_SIZE-1]; + *node = 0; + + int i; + + do + { + int child_map = tex1Dfetch(child_map_arr, *node); + int child_len = tex1Dfetch(child_len_arr, *node); + + if (*node < first_leaf) + { + for (i=0; i < child_len; i++) + if (intersect_node(origin, direction, child_map+i)) + *node++ = child_map+i; + + if (node == head) + break; + + node--; + } + else // node is a leaf + { + for (i=0; i < child_len; i++) + { + int mesh_idx = 3*(child_map + i); + + float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx)); + float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1)); + float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2)); + + if (intersect_triangle(origin, direction, v0, v1, v2, distance)) + { + if (triangle_idx == -1) + { + triangle_idx = child_map + i; + min_distance = distance; + continue; + } + + if (distance < min_distance) + { + triangle_idx = child_map + i; + min_distance = distance; + } + } + } // triangle loop + + node--; + + } // node is a leaf + + } // while loop + while (node != head); + + return triangle_idx; +} + +extern "C" +{ + +/* Translate `points` by the vector `v` */ +__global__ void translate(int max_idx, float3 *points, float3 v) +{ + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx > max_idx) + return; + + *(points+idx) += v; +} + +/* Rotate `points` through an angle `phi` counter-clockwise about the + axis `axis` (when looking towards +infinity). */ +__global__ void rotate(int max_idx, float3 *points, float phi, float3 axis) +{ + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx > max_idx) + return; + + *(points+idx) = rotate(*(points+idx), phi, axis); +} + +/* Trace the rays starting at `origins` 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 max_idx, float3 *origins, float3 *directions, int first_leaf, int *pixels) +{ + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx > max_idx) + return; + + float3 origin = *(origins+idx); + float3 direction = *(directions+idx); + direction /= norm(direction); + + int intersection_idx = intersect_mesh(origin, direction, first_leaf); + + if (intersection_idx == -1) + { + *(pixels+idx) = 0; + } + else + { + int mesh_idx = 3*intersection_idx; + + float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx)); + float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1)); + float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2)); + + *(pixels+idx) = get_color(direction, v0, v1, v2); + } +} // ray_trace + +/* Propagate the photons starting at `origins` 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 max_idx, float3 *origins, float3 *directions, int first_leaf, int *hit_solids) +{ + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx > max_idx) + return; + + float3 origin = *(origins+idx); + float3 direction = *(directions+idx); + direction /= norm(direction); + + *(hit_solids+idx) = intersect_mesh(origin, direction, first_leaf); +} // propagate + +} // extern "c" -- cgit