added a bounding volume hierarchy

2 files changed, 179 insertions, 48 deletions

diff --git a/src/intersect.cu b/src/intersect.cu
index c78350e..8c2bf45 100644
--- a/src/intersect.cu
+++ b/src/intersect.cu
@@ -1,17 +1,33 @@
 //-*-c-*-
+#include <math_constants.h>
+
+#include "linalg.h"
+#include "matrix.h"
+#include "rotate.h"
 
 texture<float4, 1, cudaReadModeElementType> mesh;
 
-__device__ bool intersect_triangle(const float3 &x, const float3 &p, const float3 &v0, const float3 &v1, const float3 &v2, float3 &intersection)
+texture<float4, 1, cudaReadModeElementType> upper_bound_arr;
+texture<float4, 1, cudaReadModeElementType> lower_bound_arr;
+
+texture<uint, 1, cudaReadModeElementType> child_map_arr;
+texture<uint, 1, cudaReadModeElementType> child_len_arr;
+
+__device__ float3 make_float3(const float4 &a)
 {
-	Matrix m = make_matrix(v1-v0, v2-v0, -p);
+	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, float3 &intersection)
+{
+	Matrix m = make_matrix(v1-v0, v2-v0, -direction);
 	
 	float determinant = det(m);
 
 	if (determinant == 0.0)
 		return false;
 
-	float3 b = x-v0;
+	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;
 
@@ -28,103 +44,213 @@ __device__ bool intersect_triangle(const float3 &x, const float3 &p, const float
 	if (u3 < 0.0 || (1-u1-u2) < 0.0)
 		return false;
 
-	intersection = x + p*u3;
+	intersection = origin + direction*u3;
 
 	return true;
 }
 
-__device__ int get_color(const float3 &p, const float3 &v0, const float3& v1, const float3 &v2)
+__device__ int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2)
 {
-	float3 normal = cross(v1-v0,v2-v0);
+	int rgb = floorf(255*dot(normalize(cross(v1-v0,v2-v0)),-direction));
+
+	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.0)
+	{
+		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.0)
+	{
+		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.0)
+	{
+		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;
+	}
 
-	float scale;
-	scale = dot(normal,-p)/(norm(normal)*norm(p));
-	if (scale < 0.0)
-		scale = dot(-normal,-p)/(norm(normal)*norm(p));
+	if (kzmax < kzmin)
+		return false;
 
-	int rgb = floorf(255*scale);
+	if (kzmin > kmin)
+		kmin = kzmin;
 
-	return rgb*65536 + rgb*256 + rgb;
+	if (kzmax < kmax)
+		kmax = kzmax;
+
+	if (kmin > kmax)
+		return false;
+
+	if (kmax < 0.0)
+		return false;
+
+	return true;
 }
 
-__device__ float3 make_float3(const float4 &a)
+__device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i)
 {
-	return make_float3(a.x, a.y, a.z);
+	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 *x, float3 v)
+__global__ void translate(int max_idx, float3 *pt, float3 v)
 {
 	int idx = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (idx > max_idx)
 		return;
 
-	x[idx] += v;
+	pt[idx] += v;
 }
 
-__global__ void rotate(int max_idx, float3 *x, float phi, float3 axis)
+__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;
 
-	x[idx] = rotate(x[idx], phi, axis);
+	pt[idx] = rotate(pt[idx], phi, axis);
 }
 
-__global__ void intersect_triangle_mesh(int max_idx, float3 *xarr, float3 *parr, int n, int *pixelarr)
+__global__ void intersect_mesh(int max_idx, float3 *origin_arr, float3 *direction_arr, int *pixel_arr, int first_leaf)
 {
 	int idx = blockIdx.x*blockDim.x + threadIdx.x;
 
 	if (idx > max_idx)
 		return;
 
-	float3 x = xarr[idx];
-	float3 p = parr[idx];
-	int *pixel = pixelarr+idx;
+	float3 origin = origin_arr[idx];
+	float3 direction = direction_arr[idx];
+
+	int *pixel = pixel_arr+idx;
 
 	bool hit = false;
 
 	float distance, min_distance;
 	float3 intersection, min_intersection;
+	
+	if (!intersect_node(origin, direction, 0))
+	{
+		*pixel = 0;
+		return;
+	}
+
+	int stack[100];
+
+	int *head = &stack[0];
+	int *node = &stack[1];
+	*node = 0;
 
 	int i;
-	for (i=0; i < n; i++)
+
+	do
 	{
-		float3 v0 = make_float3(tex1Dfetch(mesh, 3*i));
-		float3 v1 = make_float3(tex1Dfetch(mesh, 3*i+1));
-		float3 v2 = make_float3(tex1Dfetch(mesh, 3*i+2));
-			
-		if (intersect_triangle(x, p, v0, v1, v2, intersection))
-		{
-			if (!hit)
-			{
-				*pixel = get_color(p, v0, v1, v2);
+		int child_map = tex1Dfetch(child_map_arr, *node);
+		int child_len = tex1Dfetch(child_len_arr, *node);
 
-				min_distance = norm(intersection-x);
-				min_intersection = intersection;
-				hit = true;
-				continue;
-			}
+		if (*node < first_leaf)
+		{
+			for (i=0; i < child_len; i++)
+				if (intersect_node(origin, direction, child_map+i))
+					*node++ = child_map+i;
 
-			distance = norm(intersection-x);
+			if (node == head)
+				break;
+			
+			node--;
 
-			if (distance < min_distance)
+		}
+		else // node is a leaf
+		{
+			for (i=0; i < child_len; i++)
 			{
-				*pixel = get_color(p, v0, v1, v2);
+				int mesh_idx = child_map + 3*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, intersection))
+				{
+					if (!hit)
+					{
+						*pixel = get_color(direction, v0, v1, v2);
+						
+						min_distance = norm(intersection-origin);
+						min_intersection = intersection;
+						hit = true;
+						continue;
+					}
+					
+					distance = norm(intersection-origin);
+					
+					if (distance < min_distance)
+					{
+						*pixel = get_color(direction, v0, v1, v2);
+						
+						min_distance = distance;
+						min_intersection = intersection;
+					}
+				}
+
+			} // triangle loop
+
+			node--;
+
+		} // node is a leaf
+
+	} // while loop
+	while (node != head);
 
-				min_distance = distance;
-				min_intersection = intersection;
-			}
-		}
-	}
-	
 	if (!hit)
 		*pixel = 0;
-	
-}
+
+} // intersect mesh
 
 } // extern "c"
diff --git a/src/linalg.h b/src/linalg.h
index 593aa9d..fe6b1b2 100644
--- a/src/linalg.h
+++ b/src/linalg.h
@@ -113,4 +113,9 @@ __device__ __host__ float norm(const float3 &a)
 	return sqrtf(dot(a,a));
 }
 
+__device__ __host__ float3 normalize(const float3 &a)
+{
+	return a/norm(a);
+}
+
 #endif