*Added rayobject_bvh

 A bvh structure to use on the raytracer

6 files changed, 297 insertions, 48 deletions

diff --git a/source/blender/render/intern/include/rayobject.h b/source/blender/render/intern/include/rayobject.h
index d516c122bcc..3b77341f229 100644
--- a/source/blender/render/intern/include/rayobject.h
+++ b/source/blender/render/intern/include/rayobject.h
@@ -123,6 +123,12 @@ void RE_rayobject_merge_bb(RayObject *ob, float *min, float *max);
  */
 int RE_rayobject_intersect(RayObject *r, Isect *i);
 
+/*
+ * Returns distance ray must travel to hit the given bounding box
+ * BB should be in format [2][3]
+ */
+float RE_rayobject_bb_intersect(const Isect *i, const float *bb);
+
 #define ISECT_EPSILON ((float)FLT_EPSILON)
 
 
diff --git a/source/blender/render/intern/include/rayobject_rtbuild.h b/source/blender/render/intern/include/rayobject_rtbuild.h
index b1458a571dd..fbe3930149c 100644
--- a/source/blender/render/intern/include/rayobject_rtbuild.h
+++ b/source/blender/render/intern/include/rayobject_rtbuild.h
@@ -50,8 +50,8 @@ typedef struct RTBuilder
 	/* axis used (if any) on the split method */
 	int split_axis;
 	
-	/* links to child partitions calculated during splitting */
-	RayObject **child[MAX_CHILDS+1];
+	/* child partitions calculated during splitting */
+	int child_offset[MAX_CHILDS+1];
 
 } RTBuilder;
 
@@ -63,22 +63,9 @@ int rtbuild_size(RTBuilder *b);
 
 /* used during tree reorganization */
 RTBuilder* rtbuild_get_child(RTBuilder *b, int child, RTBuilder *tmp);
-void rtbuild_mean_split(RTBuilder *b, int nchilds, int axis);
-void rtbuild_mean_split_largest_axis(RTBuilder *b, int nchilds);
 
-/*
-static BVHNode *bvh_rearrange(BVHTree *tree, RTBuilder *b)
-{
-	int i;
-	int nc = rtbuild_mean_split_largest_axis(b, BVH_NCHILDS);
-	RTBuilder tmp;
-	
-	BVHNode *bvh = tree->next_node++;
-
-	bvh->split_axis = tmp->split_axis;	
-	for(i=0; i<nc; i++)
-		bvh->child[i] = bvh_rearrange( rtbuild_get_child(b, i, &tmp) );
-}
- */
+/* Calculates child partitions and returns number of efectively needed partitions */
+int rtbuild_mean_split(RTBuilder *b, int nchilds, int axis);
+int rtbuild_mean_split_largest_axis(RTBuilder *b, int nchilds);
 
 #endif
diff --git a/source/blender/render/intern/source/rayobject.c b/source/blender/render/intern/source/rayobject.c
index 3af2969b67e..7779c2aee82 100644
--- a/source/blender/render/intern/source/rayobject.c
+++ b/source/blender/render/intern/source/rayobject.c
@@ -36,6 +36,34 @@
 #include "render_types.h"
 #include "rayobject.h"
 
+
+/*
+ * Determines the distance that the ray must travel to hit the bounding volume of the given node
+ * Based on Tactical Optimization of Ray/Box Intersection, by Graham Fyffe
+ *  [http://tog.acm.org/resources/RTNews/html/rtnv21n1.html#art9]
+ */
+float RE_rayobject_bb_intersect(const Isect *isec, const float *bb)
+{
+	float dist;
+	
+	float t1x = (bb[isec->bv_index[0]] - isec->start[0]) * isec->idot_axis[0];
+	float t2x = (bb[isec->bv_index[1]] - isec->start[0]) * isec->idot_axis[0];
+	float t1y = (bb[isec->bv_index[2]] - isec->start[1]) * isec->idot_axis[1];
+	float t2y = (bb[isec->bv_index[3]] - isec->start[1]) * isec->idot_axis[1];
+	float t1z = (bb[isec->bv_index[4]] - isec->start[2]) * isec->idot_axis[2];
+	float t2z = (bb[isec->bv_index[5]] - isec->start[2]) * isec->idot_axis[2];
+
+	if(t1x > t2y || t2x < t1y || t1x > t2z || t2x < t1z || t1y > t2z || t2y < t1z) return FLT_MAX;
+	if(t2x < 0.0 || t2y < 0.0 || t2z < 0.0) return FLT_MAX;
+	if(t1x > isec->labda || t1y > isec->labda || t1z > isec->labda) return FLT_MAX;
+
+	dist = t1x;
+	if (t1y > dist) dist = t1y;
+    if (t1z > dist) dist = t1z;
+	return dist;
+}
+
+
 /* only for self-intersecting test with current render face (where ray left) */
 static int intersection2(VlakRen *face, float r0, float r1, float r2, float rx1, float ry1, float rz1)
 {
@@ -262,28 +290,43 @@ static int intersect_rayface(RayFace *face, Isect *is)
 	return 0;
 }
 
-int RE_rayobject_raycast(RayObject *r, Isect *i)
+int RE_rayobject_raycast(RayObject *r, Isect *isec)
 {
-	RE_RC_COUNT(i->count->raycast.test);
+	int i;
+	RE_RC_COUNT(isec->count->raycast.test);
+
+	/* Setup vars used on raycast */
+	isec->dist = VecLength(isec->vec);
+	
+	for(i=0; i<3; i++)
+	{
+		isec->idot_axis[i]		= 1.0f / isec->vec[i];
+		
+		isec->bv_index[2*i]		= isec->idot_axis[i] < 0.0 ? 1 : 0;
+		isec->bv_index[2*i+1]	= 1 - isec->bv_index[2*i];
+		
+		isec->bv_index[2*i]		= i+3*isec->bv_index[2*i];
+		isec->bv_index[2*i+1]	= i+3*isec->bv_index[2*i+1];
+	}
 
-	i->dist = VecLength(i->vec);
 	
-	if(i->mode==RE_RAY_SHADOW && i->last_hit && RE_rayobject_intersect(i->last_hit, i))
+	/* Last hit heuristic */
+	if(isec->mode==RE_RAY_SHADOW && isec->last_hit && RE_rayobject_intersect(isec->last_hit, isec))
 	{
-		RE_RC_COUNT(i->count->raycast.hit);
-		RE_RC_COUNT(i->count->rayshadow_last_hit_optimization );
+		RE_RC_COUNT(isec->count->raycast.hit);
+		RE_RC_COUNT(isec->count->rayshadow_last_hit_optimization );
 		return 1;
 	}
 
 #ifdef RE_RAYCOUNTER
-	if(RE_rayobject_intersect(r, i))
+	if(RE_rayobject_intersect(r, isec))
 	{
-		RE_RC_COUNT(i->count->raycast.hit);
+		RE_RC_COUNT(isec->count->raycast.hit);
 		return 1;
 	}
 	return 0;
 #else
-	return RE_rayobject_intersect(r, i);
+	return RE_rayobject_intersect(r, isec);
 #endif
 }
 
diff --git a/source/blender/render/intern/source/rayobject_bvh.c b/source/blender/render/intern/source/rayobject_bvh.c
new file mode 100644
index 00000000000..ea590ef83fd
--- /dev/null
+++ b/source/blender/render/intern/source/rayobject_bvh.c
@@ -0,0 +1,200 @@
+/**
+ * $Id$
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version. 
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2009 Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): André Pinto.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+#include <assert.h>
+
+#include "MEM_guardedalloc.h"
+#include "BKE_utildefines.h"
+#include "BLI_arithb.h"
+#include "RE_raytrace.h"
+#include "rayobject_rtbuild.h"
+#include "rayobject.h"
+
+#define BVH_NCHILDS	2
+typedef struct BVHTree BVHTree;
+
+static int  bvh_intersect(BVHTree *obj, Isect *isec);
+static void bvh_add(BVHTree *o, RayObject *ob);
+static void bvh_done(BVHTree *o);
+static void bvh_free(BVHTree *o);
+static void bvh_bb(BVHTree *o, float *min, float *max);
+
+static RayObjectAPI bvh_api =
+{
+	(RE_rayobject_raycast_callback) bvh_intersect,
+	(RE_rayobject_add_callback)     bvh_add,
+	(RE_rayobject_done_callback)    bvh_done,
+	(RE_rayobject_free_callback)    bvh_free,
+	(RE_rayobject_merge_bb_callback)bvh_bb
+};
+
+typedef struct BVHNode BVHNode;
+struct BVHNode
+{
+	BVHNode *child[BVH_NCHILDS];
+	float	 bb[2][3];
+};
+
+struct BVHTree
+{
+	RayObject rayobj;
+
+	BVHNode *alloc, *next_node, *root;
+	RTBuilder *builder;
+
+};
+
+
+RayObject *RE_rayobject_bvh_create(int size)
+{
+	BVHTree *obj= (BVHTree*)MEM_callocN(sizeof(BVHTree), "BVHTree");
+	assert( RayObject_isAligned(obj) ); /* RayObject API assumes real data to be 4-byte aligned */	
+	
+	obj->rayobj.api = &bvh_api;
+	obj->builder = rtbuild_create( size );
+	obj->root = NULL;
+	
+	return RayObject_unalignRayAPI((RayObject*) obj);
+}
+
+static void bvh_free(BVHTree *obj)
+{
+	if(obj->builder)
+		rtbuild_free(obj->builder);
+
+	if(obj->alloc)
+		MEM_freeN(obj->alloc);
+
+	MEM_freeN(obj);
+}
+
+
+static void bvh_merge_bb(BVHNode *node, float *min, float *max)
+{
+	if(RayObject_isAligned(node))
+	{
+		//TODO only half operations needed
+		DO_MINMAX(node->bb[0], min, max);
+		DO_MINMAX(node->bb[1], min, max);
+	}
+	else
+	{
+		RE_rayobject_merge_bb( (RayObject*)node, min, max);
+	}
+}
+
+static void bvh_bb(BVHTree *obj, float *min, float *max)
+{
+	bvh_merge_bb(obj->root, min, max);
+}
+
+/*
+ * Tree transverse
+ */
+static int dfs_raycast(BVHNode *node, Isect *isec)
+{
+	int hit = 0;
+	if(RE_rayobject_bb_intersect(isec, (const float*)node->bb) != FLT_MAX)
+	{
+		int i;
+		for(i=0; i<BVH_NCHILDS; i++)
+			if(RayObject_isAligned(node->child[i]))
+			{
+				hit = dfs_raycast(node->child[i], isec);
+				if(hit && isec->mode == RE_RAY_SHADOW) return hit;
+			}
+			else
+			{
+				hit = RE_rayobject_intersect( (RayObject*)node->child[i], isec);
+				if(hit && isec->mode == RE_RAY_SHADOW) return hit;
+			}
+	}
+	return hit;
+}
+
+static int bvh_intersect(BVHTree *obj, Isect *isec)
+{
+	if(RayObject_isAligned(obj->root))
+		return dfs_raycast(obj->root, isec);
+	else
+		return RE_rayobject_intersect( (RayObject*)obj->root, isec);
+}
+
+
+/*
+ * Builds a BVH tree from builder object
+ */
+static void bvh_add(BVHTree *obj, RayObject *ob)
+{
+	rtbuild_add( obj->builder, ob );
+}
+
+static BVHNode *bvh_rearrange(BVHTree *tree, RTBuilder *builder)
+{
+	if(rtbuild_size(builder) == 1)
+	{
+		return (BVHNode*)builder->begin[0];
+	}
+	else
+	{
+		int i;
+		int nc = rtbuild_mean_split_largest_axis(builder, BVH_NCHILDS);
+		RTBuilder tmp;
+	
+		BVHNode *parent = tree->next_node++;
+
+		INIT_MINMAX(parent->bb[0], parent->bb[1]);
+//		bvh->split_axis = tmp->split_axis;
+		for(i=0; i<nc; i++)
+		{
+			parent->child[i] = bvh_rearrange( tree, rtbuild_get_child(builder, i, &tmp) );
+			bvh_merge_bb(parent->child[i], parent->bb[0], parent->bb[1]);
+		}
+
+		return parent;
+	}
+}
+	
+static void bvh_done(BVHTree *obj)
+{
+	int needed_nodes;
+	assert(obj->root == NULL && obj->next_node == NULL && obj->builder);
+
+	needed_nodes = (rtbuild_size(obj->builder)+1)*2;
+	assert(needed_nodes > 0);
+	obj->alloc = (BVHNode*)MEM_mallocN( sizeof(BVHNode)*needed_nodes, "BVHTree.Nodes");
+	obj->next_node = obj->alloc;
+	obj->root = bvh_rearrange( obj, obj->builder );
+
+	assert(obj->alloc+needed_nodes >= obj->next_node);
+	
+
+	rtbuild_free( obj->builder );
+	obj->builder = NULL;
+}
+
diff --git a/source/blender/render/intern/source/rayobject_rtbuild.c b/source/blender/render/intern/source/rayobject_rtbuild.c
index b89729e2b45..46833b8c15d 100644
--- a/source/blender/render/intern/source/rayobject_rtbuild.c
+++ b/source/blender/render/intern/source/rayobject_rtbuild.c
@@ -2,12 +2,13 @@
 #include "MEM_guardedalloc.h"
 #include "BLI_arithb.h"
 #include "BKE_utildefines.h"
+#include <assert.h>
 
 static int partition_nth_element(RTBuilder *b, int _begin, int _end, int n);
 static void split_leafs(RTBuilder *b, int *nth, int partitions, int split_axis);
 
 
-static void RayObject_rtbuild_init(RTBuilder *b, RayObject **begin, RayObject **end)
+static void rtbuild_init(RTBuilder *b, RayObject **begin, RayObject **end)
 {
 	int i;
 
@@ -16,35 +17,35 @@ static void RayObject_rtbuild_init(RTBuilder *b, RayObject **begin, RayObject **
 	b->split_axis = 0;
 	
 	for(i=0; i<MAX_CHILDS; i++)
-		b->child[i] = 0;
+		b->child_offset[i] = 0;
 }
 
-RTBuilder* RayObject_rtbuild_create(int size)
+RTBuilder* rtbuild_create(int size)
 {
 	RTBuilder *builder  = (RTBuilder*) MEM_mallocN( sizeof(RTBuilder), "RTBuilder" );
-	RayObject **memblock= (RayObject**)MEM_mallocN( sizeof(RayObject*),"RTBuilder.objects");
-	RayObject_rtbuild_init(builder, memblock, memblock);
+	RayObject **memblock= (RayObject**)MEM_mallocN( sizeof(RayObject*)*size,"RTBuilder.objects");
+	rtbuild_init(builder, memblock, memblock);
 	return builder;
 }
 
-void RayObject_rtbuild_free(RTBuilder *b)
+void rtbuild_free(RTBuilder *b)
 {
 	MEM_freeN(b->begin);
 	MEM_freeN(b);
 }
 
-void RayObject_rtbuild_add(RTBuilder *b, RayObject *o)
+void rtbuild_add(RTBuilder *b, RayObject *o)
 {
 	*(b->end++) = o;
 }
 
 RTBuilder* rtbuild_get_child(RTBuilder *b, int child, RTBuilder *tmp)
 {
-	RayObject_rtbuild_init( tmp, b->child[child], b->child[child+1] );
+	rtbuild_init( tmp, b->begin + b->child_offset[child], b->begin + b->child_offset[child+1] );
 	return tmp;
 }
 
-int RayObject_rtbuild_size(RTBuilder *b)
+int rtbuild_size(RTBuilder *b)
 {
 	return b->end - b->begin;
 }
@@ -86,17 +87,23 @@ static int calc_largest_axis(RTBuilder *b)
 //Unballanced mean
 //TODO better balance nodes
 //TODO suport for variable number of partitions (its hardcoded in 2)
-void rtbuild_mean_split(RTBuilder *b, int nchilds, int axis)
+int rtbuild_mean_split(RTBuilder *b, int nchilds, int axis)
 {
-	int nth[3] = {0, (b->end - b->begin)/2, b->end-b->begin};
-	split_leafs(b, nth, 2, axis);
+	b->child_offset[0] = 0;
+	b->child_offset[1] = (b->end - b->begin) / 2;
+	b->child_offset[2] = (b->end - b->begin);
+	
+	assert( b->child_offset[0] != b->child_offset[1] && b->child_offset[1] != b->child_offset[2]);
+
+	split_leafs(b, b->child_offset, 2, axis);
+	return 2;
 }
 	
 	
-void rtbuild_mean_split_largest_axis(RTBuilder *b, int nchilds)
+int rtbuild_mean_split_largest_axis(RTBuilder *b, int nchilds)
 {
 	int axis = calc_largest_axis(b);
-	rtbuild_mean_split(b, nchilds, axis);
+	return rtbuild_mean_split(b, nchilds, axis);
 }
 
 
@@ -110,11 +117,12 @@ static void sort_swap(RTBuilder *b, int i, int j)
 	SWAP(RayObject*, b->begin[i], b->begin[j]);
 }
  
-static int sort_get_value(RTBuilder *b, int i)
+static float sort_get_value(RTBuilder *b, int i)
 {
 	float min[3], max[3];
+	INIT_MINMAX(min, max);
 	RE_rayobject_merge_bb(b->begin[i], min, max);
-	return max[i];
+	return max[b->split_axis];
 }
  
 static int medianof3(RTBuilder *d, int a, int b, int c)
@@ -174,9 +182,14 @@ static int partition(RTBuilder *b, int lo, int mid, int hi)
 	int i=lo, j=hi;
 	while (1)
 	{
-		while (sort_get_value(b,i) < x) i++;
+		while(sort_get_value(b,i) < x)
+			i++;
+
 		j--;
-		while (x < sort_get_value(b,j)) j--;
+
+		while(x < sort_get_value(b,j))
+			j--;
+
 		if(!(i < j))
 			return i;
 
diff --git a/source/blender/render/intern/source/rayshade.c b/source/blender/render/intern/source/rayshade.c
index 75c93fc1c80..4c16ac05ec2 100644
--- a/source/blender/render/intern/source/rayshade.c
+++ b/source/blender/render/intern/source/rayshade.c
@@ -180,7 +180,7 @@ RayObject* makeraytree_object(Render *re, ObjectInstanceRen *obi)
 		if(re->r.raystructure == R_RAYSTRUCTURE_HIER_BVH_OCTREE)
 			raytree = obr->raytree = RE_rayobject_octree_create( re->r.ocres, faces );
 		else //if(re->r.raystructure == R_RAYSTRUCTURE_HIER_BVH_BVH)
-			raytree = obr->raytree = RE_rayobject_blibvh_create( faces );
+			raytree = obr->raytree = RE_rayobject_tree_create( faces );
 			
 		face = obr->rayfaces = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "ObjectRen faces");
 		obr->rayobi = obi;
@@ -240,7 +240,7 @@ static void makeraytree_hier(Render *re)
 		num_objects++;
 
 	//Create raytree
-	re->raytree = RE_rayobject_blibvh_create( num_objects );
+	re->raytree = RE_rayobject_tree_create( num_objects );
 	
 	for(obi=re->instancetable.first; obi; obi=obi->next)
 	if(is_raytraceable(re, obi))
@@ -292,7 +292,7 @@ static void makeraytree_single(Render *re)
 	if(re->r.raystructure == R_RAYSTRUCTURE_SINGLE_OCTREE)
 		raytree = re->raytree = RE_rayobject_octree_create( re->r.ocres, faces );
 	else //if(re->r.raystructure == R_RAYSTRUCTURE_SINGLE_BVH)
-		raytree = re->raytree = RE_rayobject_blibvh_create( faces );
+		raytree = re->raytree = RE_rayobject_tree_create( faces );
 
 	face	= re->rayfaces	= (RayFace*)MEM_callocN(faces*sizeof(RayFace), "Render ray faces");