*Moved rtbuild to bf_render_raytrace

*Added vbvh - Just a experimental tree type :)
Variable Way BVH - there is no hardcoded number of childs per each Tree Node
 - idea is to optimize a tree to reduced the expected number of BB tests even after applying SAH (for that an hardcoded n-way is not enough)
 - for now childs are stored on a linked list

4 files changed, 789 insertions, 8 deletions

diff --git a/source/blender/render/intern/raytrace/bvh.h b/source/blender/render/intern/raytrace/bvh.h
index 44f531faa9f..d9277f9a583 100644
--- a/source/blender/render/intern/raytrace/bvh.h
+++ b/source/blender/render/intern/raytrace/bvh.h
@@ -99,7 +99,12 @@ static int bvh_node_stack_raycast(Node *root, Isect *isec)
 	Node *stack[MAX_STACK_SIZE];
 	int hit = 0, stack_pos = 0;
 		
-	stack[stack_pos++] = root;
+	//Assume the BB of root always succeed
+	if(1)
+		bvh_node_push_childs(root, isec, stack, stack_pos);
+	else
+		stack[stack_pos++] = root;
+
 	while(stack_pos)
 	{
 		Node *node = stack[--stack_pos];
diff --git a/source/blender/render/intern/raytrace/rayobject_bvh.cpp b/source/blender/render/intern/raytrace/rayobject_bvh.cpp
index 8a63e088a34..2c2a260df98 100644
--- a/source/blender/render/intern/raytrace/rayobject_bvh.cpp
+++ b/source/blender/render/intern/raytrace/rayobject_bvh.cpp
@@ -26,18 +26,15 @@
  *
  * ***** END GPL LICENSE BLOCK *****
  */
-extern "C"
-{
 #include <assert.h>
+
+#include "RE_raytrace.h"
+#include "rayobject_rtbuild.h"
+#include "rayobject.h"
 #include "MEM_guardedalloc.h"
 #include "BKE_utildefines.h"
 #include "BLI_arithb.h"
 #include "BLI_memarena.h"
-#include "RE_raytrace.h"
-#include "rayobject_rtbuild.h"
-#include "rayobject.h"
-};
-
 #include "bvh.h"
 
 #define BVH_NCHILDS 2
diff --git a/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp b/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp
new file mode 100644
index 00000000000..dff0b02c00e
--- /dev/null
+++ b/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp
@@ -0,0 +1,531 @@
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "rayobject_rtbuild.h"
+#include "MEM_guardedalloc.h"
+#include "BLI_arithb.h"
+#include "BKE_utildefines.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 int split_leafs_by_plane(RTBuilder *b, int begin, int end, float plane);
+
+static void rtbuild_init(RTBuilder *b, RayObject **begin, RayObject **end)
+{
+	int i;
+
+	b->begin = begin;
+	b->end   = end;
+	b->split_axis = 0;
+	b->child_sorted_axis = -1;
+	
+	for(i=0; i<RTBUILD_MAX_CHILDS; i++)
+		b->child_offset[i] = 0;
+		
+	INIT_MINMAX(b->bb, b->bb+3);
+}
+
+RTBuilder* rtbuild_create(int size)
+{
+	RTBuilder *builder  = (RTBuilder*) MEM_mallocN( sizeof(RTBuilder), "RTBuilder" );
+	RayObject **memblock= (RayObject**)MEM_mallocN( sizeof(RayObject*)*size,"RTBuilder.objects");
+	rtbuild_init(builder, memblock, memblock);
+	return builder;
+}
+
+void rtbuild_free(RTBuilder *b)
+{
+	MEM_freeN(b->begin);
+	MEM_freeN(b);
+}
+
+void rtbuild_add(RTBuilder *b, RayObject *o)
+{
+	*(b->end++) = o;
+}
+
+void rtbuild_calc_bb(RTBuilder *b)
+{
+	if(b->bb[0] == 1.0e30f)
+	{
+		RayObject **index = b->begin;
+		for(; index != b->end; index++)
+			RE_rayobject_merge_bb(*index, b->bb, b->bb+3);
+	}
+}
+
+void rtbuild_merge_bb(RTBuilder *b, float *min, float *max)
+{
+	rtbuild_calc_bb(b);
+	DO_MIN(b->bb, min);
+	DO_MAX(b->bb+3, max);
+}
+
+int rtbuild_get_largest_axis(RTBuilder *b)
+{
+	rtbuild_calc_bb(b);
+	return bb_largest_axis(b->bb, b->bb+3);
+}
+
+
+int rtbuild_size(RTBuilder *b)
+{
+	return b->end - b->begin;
+}
+
+
+RTBuilder* rtbuild_get_child(RTBuilder *b, int child, RTBuilder *tmp)
+{
+	rtbuild_init( tmp, b->begin + b->child_offset[child], b->begin + b->child_offset[child+1] );
+	tmp->child_sorted_axis = b->child_sorted_axis;
+	return tmp;
+}
+
+
+
+//Left balanced tree
+int rtbuild_mean_split(RTBuilder *b, int nchilds, int axis)
+{
+	int i;
+	int mleafs_per_child, Mleafs_per_child;
+	int tot_leafs  = rtbuild_size(b);
+	int missing_leafs;
+
+	long long s;
+
+	assert(nchilds <= RTBUILD_MAX_CHILDS);
+	
+	//TODO optimize calc of leafs_per_child
+	for(s=nchilds; s<tot_leafs; s*=nchilds);
+	Mleafs_per_child = s/nchilds;
+	mleafs_per_child = Mleafs_per_child/nchilds;
+	
+	//split min leafs per child	
+	b->child_offset[0] = 0;
+	for(i=1; i<=nchilds; i++)
+		b->child_offset[i] = mleafs_per_child;
+	
+	//split remaining leafs
+	missing_leafs = tot_leafs - mleafs_per_child*nchilds;
+	for(i=1; i<=nchilds; i++)
+	{
+		if(missing_leafs > Mleafs_per_child - mleafs_per_child)
+		{
+			b->child_offset[i] += Mleafs_per_child - mleafs_per_child;
+			missing_leafs -= Mleafs_per_child - mleafs_per_child;
+		}
+		else
+		{
+			b->child_offset[i] += missing_leafs;
+			missing_leafs = 0;
+			break;
+		}
+	}
+	
+	//adjust for accumulative offsets
+	for(i=1; i<=nchilds; i++)
+		b->child_offset[i] += b->child_offset[i-1];
+
+	//Count created childs
+	for(i=nchilds; b->child_offset[i] == b->child_offset[i-1]; i--);
+	split_leafs(b, b->child_offset, i, axis);
+	
+	assert( b->child_offset[0] == 0 && b->child_offset[i] == tot_leafs );
+	return i;
+}
+	
+	
+int rtbuild_mean_split_largest_axis(RTBuilder *b, int nchilds)
+{
+	int axis = rtbuild_get_largest_axis(b);
+	return rtbuild_mean_split(b, nchilds, axis);
+}
+
+/*
+ * "separators" is an array of dim NCHILDS-1
+ * and indicates where to cut the childs
+ */
+int rtbuild_median_split(RTBuilder *b, float *separators, int nchilds, int axis)
+{
+	int size = rtbuild_size(b);
+		
+	assert(nchilds <= RTBUILD_MAX_CHILDS);
+	if(size <= nchilds)
+	{
+		return rtbuild_mean_split(b, nchilds, axis);
+	}
+	else
+	{
+		int i;
+
+		b->split_axis = axis;
+		
+		//Calculate child offsets
+		b->child_offset[0] = 0;
+		for(i=0; i<nchilds-1; i++)
+			b->child_offset[i+1] = split_leafs_by_plane(b, b->child_offset[i], size, separators[i]);
+		b->child_offset[nchilds] = size;
+		
+		for(i=0; i<nchilds; i++)
+			if(b->child_offset[i+1] - b->child_offset[i] == size)
+				return rtbuild_mean_split(b, nchilds, axis);
+		
+		return nchilds;
+	}	
+}
+
+int rtbuild_median_split_largest_axis(RTBuilder *b, int nchilds)
+{
+	int la, i;
+	float separators[RTBUILD_MAX_CHILDS];
+	
+	rtbuild_calc_bb(b);
+
+	la = bb_largest_axis(b->bb,b->bb+3);
+	for(i=1; i<nchilds; i++)
+		separators[i-1] = (b->bb[la+3]-b->bb[la])*i / nchilds;
+		
+	return rtbuild_median_split(b, separators, nchilds, la);
+}
+
+//Heuristics Object Splitter
+typedef struct CostObject CostObject;
+struct CostObject
+{
+	RayObject *obj;
+	float cost;
+	float bb[6];
+};
+
+//Ugly.. but using qsort and no globals its the cleaner I can get
+#define costobject_cmp(axis) costobject_cmp##axis
+#define define_costobject_cmp(axis) 								\
+int costobject_cmp(axis)(const CostObject *a, const CostObject *b)	\
+{																	\
+	if(a->bb[axis] < b->bb[axis]) return -1;						\
+	if(a->bb[axis] > b->bb[axis]) return  1;						\
+	if(a->obj < b->obj) return -1;									\
+	if(a->obj > b->obj) return 1;									\
+	return 0;														\
+}	
+
+define_costobject_cmp(0)
+define_costobject_cmp(1)
+define_costobject_cmp(2)
+define_costobject_cmp(3)
+define_costobject_cmp(4)
+define_costobject_cmp(5)
+
+void costobject_sort(CostObject *begin, CostObject *end, int axis)
+{
+	//TODO introsort
+	     if(axis == 0) qsort(begin, end-begin, sizeof(*begin), (int(*)(const void *, const void *)) costobject_cmp(0));
+	else if(axis == 1) qsort(begin, end-begin, sizeof(*begin), (int(*)(const void *, const void *)) costobject_cmp(1));
+	else if(axis == 2) qsort(begin, end-begin, sizeof(*begin), (int(*)(const void *, const void *)) costobject_cmp(2));
+	else if(axis == 3) qsort(begin, end-begin, sizeof(*begin), (int(*)(const void *, const void *)) costobject_cmp(3));
+	else if(axis == 4) qsort(begin, end-begin, sizeof(*begin), (int(*)(const void *, const void *)) costobject_cmp(4));
+	else if(axis == 5) qsort(begin, end-begin, sizeof(*begin), (int(*)(const void *, const void *)) costobject_cmp(5));
+}
+
+
+
+/* Object Surface Area Heuristic splitter */
+int rtbuild_heuristic_object_split(RTBuilder *b, int nchilds)
+{
+	int size = rtbuild_size(b);		
+	assert(nchilds == 2);
+
+	if(size <= nchilds)
+	{
+		return rtbuild_mean_split_largest_axis(b, nchilds);
+	}
+	else
+	{
+		float bcost = FLT_MAX;
+		float childrens_cost = 0;
+		int i, axis, baxis = -1, boffset = size/2, k, try_axis[3];
+		CostObject *cost   = (CostObject*)MEM_mallocN( sizeof(CostObject)*size, "RTBuilder.HeuristicObjectSplitter" );
+		float      *acc_bb = (float*)MEM_mallocN( sizeof(float)*6*size, "RTBuilder.HeuristicObjectSplitterBB" );
+
+		for(i=0; i<size; i++)
+		{
+			cost[i].obj = b->begin[i];
+			INIT_MINMAX(cost[i].bb, cost[i].bb+3);
+			RE_rayobject_merge_bb(cost[i].obj, cost[i].bb, cost[i].bb+3);
+			cost[i].cost = RE_rayobject_cost(cost[i].obj);
+			childrens_cost += cost[i].cost;
+		}
+		
+		if(b->child_sorted_axis >= 0 && b->child_sorted_axis < 3)
+		{
+			try_axis[0] =  b->child_sorted_axis;
+			try_axis[1] = (b->child_sorted_axis+1)%3;
+			try_axis[2] = (b->child_sorted_axis+2)%3;
+		}
+		else
+		{
+			try_axis[0] = 0;
+			try_axis[1] = 1;
+			try_axis[2] = 2;
+		}
+		
+		for(axis=try_axis[k=0]; k<3; axis=try_axis[++k])
+		{
+			float left_cost, right_cost;
+			float other_bb[6];
+
+
+
+			costobject_sort(cost, cost+size, axis);
+			for(i=size-1; i>=0; i--)
+			{
+				float *bb = acc_bb+i*6;
+				if(i == size-1)
+				{
+					VECCOPY(bb, cost[i].bb);
+					VECCOPY(bb+3, cost[i].bb+3);
+				}
+				else
+				{
+					bb[0] = MIN2(cost[i].bb[0], bb[6+0]);
+					bb[1] = MIN2(cost[i].bb[1], bb[6+1]);
+					bb[2] = MIN2(cost[i].bb[2], bb[6+2]);
+
+					bb[3] = MAX2(cost[i].bb[3], bb[6+3]);
+					bb[4] = MAX2(cost[i].bb[4], bb[6+4]);
+					bb[5] = MAX2(cost[i].bb[5], bb[6+5]);
+				}
+			}
+			
+			INIT_MINMAX(other_bb, other_bb+3);
+			DO_MIN( cost[0].bb,   other_bb   );
+			DO_MAX( cost[0].bb+3, other_bb+3 );
+			left_cost  = cost[0].cost;
+			right_cost = childrens_cost-cost[0].cost;
+			if(right_cost < 0) right_cost = 0;
+
+			for(i=1; i<size; i++)
+			{
+				//Worst case heuristic (cost of each child is linear)
+				float hcost, left_side, right_side;
+				
+				left_side = bb_area(other_bb, other_bb+3)*(left_cost+logf(i));
+				right_side= bb_area(acc_bb+i*6, acc_bb+i*6+3)*(right_cost+logf(i));
+				
+				if(left_side > bcost) break;	//No way we can find a better heuristic in this axis
+
+				hcost = left_side+right_side;
+				assert(hcost >= 0);
+				if( hcost < bcost
+				|| (hcost == bcost && axis < baxis)) //this makes sure the tree built is the same whatever is the order of the sorting axis
+				{
+					bcost = hcost;
+					baxis = axis;
+					boffset = i;
+				}
+				DO_MIN( cost[i].bb,   other_bb   );
+				DO_MAX( cost[i].bb+3, other_bb+3 );
+				left_cost  += cost[i].cost;
+				right_cost -= cost[i].cost;
+				if(right_cost < 0.0f) right_cost = 0.0;
+			}
+			
+			if(baxis == axis)
+			{
+				for(i=0; i<size; i++)
+					b->begin[i] = cost[i].obj;
+				b->child_sorted_axis = axis;
+			}
+			
+			assert(baxis >= 0 && baxis < 3);
+		}
+			
+		b->child_offset[0] = 0;
+		b->child_offset[1] = boffset;
+		b->child_offset[2] = size;
+		
+		MEM_freeN(acc_bb);
+		MEM_freeN(cost);
+		return nchilds;		
+	}
+}
+
+/*
+ * Helper code
+ * PARTITION code / used on mean-split
+ * basicly this a std::nth_element (like on C++ STL algorithm)
+ */
+static void sort_swap(RTBuilder *b, int i, int j)
+{
+	SWAP(RayObject*, b->begin[i], b->begin[j]);
+}
+ 
+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[b->split_axis];
+}
+ 
+static int medianof3(RTBuilder *d, int a, int b, int c)
+{
+	float fa = sort_get_value( d, a );
+	float fb = sort_get_value( d, b );
+	float fc = sort_get_value( d, c );
+
+	if(fb < fa)
+	{
+		if(fc < fb)
+			return b;
+		else
+		{
+			if(fc < fa)
+				return c;
+			else
+				return a;
+		}
+	}
+	else
+	{
+		if(fc < fb)
+		{
+			if(fc < fa)
+				return a;
+			else
+				return c;
+		}
+		else
+			return b;
+	}
+}
+
+static void insertionsort(RTBuilder *b, int lo, int hi)
+{
+	int i;
+	for(i=lo; i<hi; i++)
+	{
+		RayObject *t = b->begin[i];
+		float tv= sort_get_value(b, i);
+		int j=i;
+		
+		while( j != lo && tv < sort_get_value(b, j-1))
+		{
+			b->begin[j] = b->begin[j-1];
+			j--;
+		}
+		b->begin[j] = t;
+	}
+}
+
+static int partition(RTBuilder *b, int lo, int mid, int hi)
+{
+	float x = sort_get_value( b, mid );
+	
+	int i=lo, j=hi;
+	while (1)
+	{
+		while(sort_get_value(b,i) < x) i++;
+		j--;
+		while(x < sort_get_value(b,j)) j--;
+		if(!(i < j))
+			return i;
+
+		sort_swap(b, i, j);
+		i++;
+	}
+}
+
+//
+// PARTITION code / used on mean-split
+// basicly this is an adapted std::nth_element (C++ STL <algorithm>)
+//
+// after a call to this function you can expect one of:
+//      every node to left of a[n] are smaller or equal to it
+//      every node to the right of a[n] are greater or equal to it
+static int partition_nth_element(RTBuilder *b, int _begin, int n, int _end)
+{
+	int begin = _begin, end = _end, cut;
+	while(end-begin > 3)
+	{
+		cut = partition(b, begin, medianof3(b, begin, begin+(end-begin)/2, end-1), end);
+		if(cut <= n)
+			begin = cut;
+		else
+			end = cut;
+	}
+	insertionsort(b, begin, end);
+
+	return n;
+}
+
+static void split_leafs(RTBuilder *b, int *nth, int partitions, int split_axis)
+{
+	int i;
+	b->split_axis = split_axis;
+
+	for(i=0; i < partitions-1; i++)
+	{
+		assert(nth[i] < nth[i+1] && nth[i+1] < nth[partitions]);
+
+		partition_nth_element(b, nth[i],  nth[i+1], nth[partitions] );
+	}
+}
+
+static int split_leafs_by_plane(RTBuilder *b, int begin, int end, float plane)
+{
+	int i;
+	for(i = begin; i != end; i++)
+	{
+		if(sort_get_value(b, i) < plane)
+		{
+			sort_swap(b, i, begin);
+			begin++;
+		}
+	}
+	return begin;
+}
+
+/*
+ * Bounding Box utils
+ */
+float bb_volume(float *min, float *max)
+{
+	return (max[0]-min[0])*(max[1]-min[1])*(max[2]-min[2]);
+}
+
+float bb_area(float *min, float *max)
+{
+	float sub[3], a;
+	sub[0] = max[0]-min[0];
+	sub[1] = max[1]-min[1];
+	sub[2] = max[2]-min[2];
+
+	a = (sub[0]*sub[1] + sub[0]*sub[2] + sub[1]*sub[2])*2;
+	assert(a >= 0.0);
+	return a;
+}
+
+int bb_largest_axis(float *min, float *max)
+{
+	float sub[3];
+	
+	sub[0] = max[0]-min[0];
+	sub[1] = max[1]-min[1];
+	sub[2] = max[2]-min[2];
+	if(sub[0] > sub[1])
+	{
+		if(sub[0] > sub[2])
+			return 0;
+		else
+			return 2;
+	}
+	else
+	{
+		if(sub[1] > sub[2])
+			return 1;
+		else
+			return 2;
+	}	
+}
diff --git a/source/blender/render/intern/raytrace/rayobject_vbvh.cpp b/source/blender/render/intern/raytrace/rayobject_vbvh.cpp
new file mode 100644
index 00000000000..3f8f69a5abe
--- /dev/null
+++ b/source/blender/render/intern/raytrace/rayobject_vbvh.cpp
@@ -0,0 +1,248 @@
+/**
+ * $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 *****
+ */
+extern "C"
+{
+#include <assert.h>
+#include "MEM_guardedalloc.h"
+#include "BKE_utildefines.h"
+#include "BLI_arithb.h"
+#include "BLI_memarena.h"
+#include "RE_raytrace.h"
+#include "rayobject_rtbuild.h"
+#include "rayobject.h"
+};
+
+#include "bvh.h"
+#include <queue>
+
+#define BVHNode VBVHNode
+#define BVHTree VBVHTree
+
+#define RAY_BB_TEST_COST (0.2f)
+#define DFS_STACK_SIZE	128
+#define DYNAMIC_ALLOC
+
+//#define rtbuild_split	rtbuild_mean_split_largest_axis		/* objects mean split on the longest axis, childs BB are allowed to overlap */
+//#define rtbuild_split	rtbuild_median_split_largest_axis	/* space median split on the longest axis, childs BB are allowed to overlap */
+#define rtbuild_split	rtbuild_heuristic_object_split		/* split objects using heuristic */
+
+struct BVHNode
+{
+	BVHNode *child;
+	BVHNode *sibling;
+
+	float	bb[6];
+};
+
+struct BVHTree
+{
+	RayObject rayobj;
+
+	BVHNode *root;
+
+	MemArena *node_arena;
+
+	float cost;
+	RTBuilder *builder;
+};
+
+
+/*
+ * Push nodes (used on dfs)
+ */
+template<class Node>
+inline static void bvh_node_push_childs(Node *node, Isect *isec, Node **stack, int &stack_pos)
+{
+	Node *child = node->child;
+	while(child)
+	{
+		stack[stack_pos++] = child;
+		if(RayObject_isAligned(child))
+			child = child->sibling;
+		else break;
+	}
+}
+
+/*
+ * BVH done
+ */
+static BVHNode *bvh_new_node(BVHTree *tree)
+{
+	BVHNode *node = (BVHNode*)BLI_memarena_alloc(tree->node_arena, sizeof(BVHNode));
+	node->sibling = NULL;
+	node->child   = NULL;
+
+	assert(RayObject_isAligned(node));
+	return node;
+}
+
+template<class Builder>
+float rtbuild_area(Builder *builder)
+{
+	float min[3], max[3];
+	INIT_MINMAX(min, max);
+	rtbuild_merge_bb(builder, min, max);
+	return bb_area(min, max);	
+}
+
+template<class Node>
+void bvh_update_bb(Node *node)
+{
+	INIT_MINMAX(node->bb, node->bb+3);
+	Node *child = node->child;
+	
+	while(child)
+	{
+		bvh_node_merge_bb(child, node->bb, node->bb+3);
+		if(RayObject_isAligned(child))
+			child = child->sibling;
+		else
+			child = 0;
+	}
+}
+
+
+template<class Tree, class Node, class Builder>
+Node *bvh_rearrange(Tree *tree, Builder *builder, float *cost)
+{
+	
+	int size = rtbuild_size(builder);
+	if(size == 1)
+	{
+		Node *node = bvh_new_node(tree);
+		INIT_MINMAX(node->bb, node->bb+3);
+		rtbuild_merge_bb(builder, node->bb, node->bb+3);
+		
+		node->child = (BVHNode*)builder->begin[0];
+
+		*cost = RE_rayobject_cost((RayObject*)node->child)+RAY_BB_TEST_COST;
+		return node;
+	}
+	else
+	{
+		Node *node = bvh_new_node(tree);
+		float parent_area;
+		
+		INIT_MINMAX(node->bb, node->bb+3);
+		rtbuild_merge_bb(builder, node->bb, node->bb+3);
+		
+		parent_area = bb_area( node->bb, node->bb+3 );
+		Node **child = &node->child;
+		
+		std::queue<Builder> childs;
+		childs.push(*builder);
+		
+		*cost = 0;
+		
+		while(!childs.empty())
+		{
+			Builder b = childs.front();
+						childs.pop();
+			
+			float hit_prob = rtbuild_area(&b) / parent_area;
+			if(hit_prob > 1.0f / 2.0f && rtbuild_size(&b) > 1)
+			{
+				//The expected number of BB test is smaller if we directly add the 2 childs of this node
+				int nc = rtbuild_split(&b, 2);
+				assert(nc == 2);
+				for(int i=0; i<nc; i++)
+				{
+					Builder tmp;
+					rtbuild_get_child(&b, i, &tmp);
+					childs.push(tmp);
+				}
+				
+			}
+			else
+			{
+				float tcost;
+				*child = bvh_rearrange<Tree,Node,Builder>(tree, &b, &tcost);
+				child = &((*child)->sibling);
+				
+				*cost += tcost*hit_prob + RAY_BB_TEST_COST;
+			}
+		}
+		assert(child != &node->child);
+		*child = 0;
+
+		return node;
+	}
+}
+
+template<>
+void bvh_done<BVHTree>(BVHTree *obj)
+{
+	int needed_nodes = (rtbuild_size(obj->builder)+1)*2;
+	if(needed_nodes > BLI_MEMARENA_STD_BUFSIZE)
+		needed_nodes = BLI_MEMARENA_STD_BUFSIZE;
+
+	obj->node_arena = BLI_memarena_new(needed_nodes);
+	BLI_memarena_use_malloc(obj->node_arena);
+
+	
+	obj->root = bvh_rearrange<BVHTree,BVHNode,RTBuilder>( obj, obj->builder, &obj->cost );
+	obj->cost = 1.0;
+	
+	rtbuild_free( obj->builder );
+	obj->builder = NULL;
+}
+
+template<>
+int bvh_intersect<BVHTree>(BVHTree *obj, Isect* isec)
+{
+	if(RayObject_isAligned(obj->root))
+		return bvh_node_stack_raycast<BVHNode,DFS_STACK_SIZE>(obj->root, isec);
+	else
+		return RE_rayobject_intersect( (RayObject*) obj->root, isec );
+}
+
+/* the cast to pointer function is needed to workarround gcc bug: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=11407 */
+static RayObjectAPI bvh_api =
+{
+	(RE_rayobject_raycast_callback) ((int(*)(BVHTree*,Isect*)) &bvh_intersect<BVHTree>),
+	(RE_rayobject_add_callback)     ((void(*)(BVHTree*,RayObject*)) &bvh_add<BVHTree>),
+	(RE_rayobject_done_callback)    ((void(*)(BVHTree*))       &bvh_done<BVHTree>),
+	(RE_rayobject_free_callback)    ((void(*)(BVHTree*))       &bvh_free<BVHTree>),
+	(RE_rayobject_merge_bb_callback)((void(*)(BVHTree*,float*,float*)) &bvh_bb<BVHTree>),
+	(RE_rayobject_cost_callback)	((float(*)(BVHTree*))      &bvh_cost<BVHTree>)
+};
+
+RayObject *RE_rayobject_vbvh_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->root = NULL;
+	
+	obj->node_arena = NULL;
+	obj->builder    = rtbuild_create( size );
+	
+	return RayObject_unalignRayAPI((RayObject*) obj);
+}