From 909309302c83423994e9c1dd36a3309890a67b90 Mon Sep 17 00:00:00 2001 From: Anthony LaTorre Date: Tue, 17 May 2011 16:58:36 -0400 Subject: added documentation --- src/intersect.cu | 300 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/kernel.cu | 270 ------------------------------------------------- 2 files changed, 300 insertions(+), 270 deletions(-) create mode 100644 src/intersect.cu delete mode 100644 src/kernel.cu (limited to 'src') diff --git a/src/intersect.cu b/src/intersect.cu new file mode 100644 index 0000000..1402aa1 --- /dev/null +++ b/src/intersect.cu @@ -0,0 +1,300 @@ +//-*-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/kernel.cu b/src/kernel.cu deleted file mode 100644 index c2b3fb2..0000000 --- a/src/kernel.cu +++ /dev/null @@ -1,270 +0,0 @@ -//-*-c-*- -#include - -#include "linalg.h" -#include "matrix.h" -#include "rotate.h" - -texture mesh; - -texture upper_bound_arr; -texture lower_bound_arr; - -texture child_map_arr; -texture child_len_arr; - -__device__ float3 make_float3(const float4 &a) -{ - return make_float3(a.x, a.y, a.z); -} - -__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 = norm(direction*u3); - - return true; -} - -__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; -} - - -__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; -} - -__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 = mesh_idx; - - min_distance = distance; - - hit = true; - continue; - } - - if (distance < min_distance) - { - *pixel = get_color(direction, v0, v1, v2); - *state = mesh_idx; - - min_distance = distance; - } - } - - } // triangle loop - - node--; - - } // node is a leaf - - } // while loop - while (node != head); - - if (!hit) - *pixel = 0; - -} // intersect mesh - -} // extern "c" -- cgit