NlogN building:

	sort once
	select subsets and kept the order (on X, Y and Z)

1 files changed, 182 insertions, 247 deletions

diff --git a/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp b/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp
index b0ac54bef2c..fb2d05a0a14 100644
--- a/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp
+++ b/source/blender/render/intern/raytrace/rayobject_rtbuild.cpp
@@ -8,21 +8,23 @@
 #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)
+static bool selected_node(RTBuilder::Object *node)
 {
-	int i;
+	return node->selected;
+}
 
-	b->begin = begin;
-	b->end   = end;
-	b->split_axis = 0;
-	b->child_sorted_axis = -1;
+static void rtbuild_init(RTBuilder *b)
+{
+	b->split_axis = -1;
+	b->primitives.begin   = 0;
+	b->primitives.end     = 0;
+	b->primitives.maxsize = 0;
 	
-	for(i=0; i<RTBUILD_MAX_CHILDS; i++)
+	for(int i=0; i<RTBUILD_MAX_CHILDS; i++)
 		b->child_offset[i] = 0;
+
+	for(int i=0; i<3; i++)
+		b->sorted_begin[i] = b->sorted_end[i] = 0;
 		
 	INIT_MINMAX(b->bb, b->bb+3);
 }
@@ -30,60 +32,129 @@ static void rtbuild_init(RTBuilder *b, RayObject **begin, RayObject **end)
 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);
+	RTBuilder::Object *memblock= (RTBuilder::Object*)MEM_mallocN( sizeof(RTBuilder::Object)*size,"RTBuilder.objects");
+
+
+	rtbuild_init(builder);
+	
+	builder->primitives.begin = builder->primitives.end = memblock;
+	builder->primitives.maxsize = size;
+	
+	for(int i=0; i<3; i++)
+	{
+		builder->sorted_begin[i] = (RTBuilder::Object**)MEM_mallocN( sizeof(RTBuilder::Object*)*size,"RTBuilder.sorted_objects");
+		builder->sorted_end[i]   = builder->sorted_begin[i];
+	} 
+	
+
 	return builder;
 }
 
 void rtbuild_free(RTBuilder *b)
 {
-	MEM_freeN(b->begin);
+	if(b->primitives.begin) MEM_freeN(b->primitives.begin);
+
+	for(int i=0; i<3; i++)
+		if(b->sorted_begin[i])
+			MEM_freeN(b->sorted_begin[i]);
+
 	MEM_freeN(b);
 }
 
 void rtbuild_add(RTBuilder *b, RayObject *o)
 {
-	*(b->end++) = o;
-}
+	assert( b->primitives.begin + b->primitives.maxsize != b->primitives.end );
+	
+	b->primitives.end->obj = o;
+	b->primitives.end->cost = RE_rayobject_cost(o);
 
-void rtbuild_calc_bb(RTBuilder *b)
-{
-	if(b->bb[0] == 1.0e30f)
+	INIT_MINMAX(b->primitives.end->bb, b->primitives.end->bb+3);
+	RE_rayobject_merge_bb(o, b->primitives.end->bb, b->primitives.end->bb+3);
+	
+	for(int i=0; i<3; i++)
 	{
-		RayObject **index = b->begin;
-		for(; index != b->end; index++)
-			RE_rayobject_merge_bb(*index, b->bb, b->bb+3);
+		*(b->sorted_end[i]) = b->primitives.end;
+		b->sorted_end[i]++;
 	}
+	b->primitives.end++;
 }
 
-void rtbuild_merge_bb(RTBuilder *b, float *min, float *max)
+int rtbuild_size(RTBuilder *b)
 {
-	rtbuild_calc_bb(b);
-	DO_MIN(b->bb, min);
-	DO_MAX(b->bb+3, max);
+	return b->sorted_end[0] - b->sorted_begin[0];
 }
 
-int rtbuild_get_largest_axis(RTBuilder *b)
+
+template<class Obj,int Axis>
+static bool obj_bb_compare(const Obj &a, const Obj &b)
 {
-	rtbuild_calc_bb(b);
-	return bb_largest_axis(b->bb, b->bb+3);
+	if(a->bb[Axis] != b->bb[Axis])
+		return a->bb[Axis] < b->bb[Axis];
+	return a->obj < b->obj;
 }
 
+template<class Item>
+static void object_sort(Item *begin, Item *end, int axis)
+{
+	if(axis == 0) return std::sort(begin, end, obj_bb_compare<Item,0> );
+	if(axis == 1) return std::sort(begin, end, obj_bb_compare<Item,1> );
+	if(axis == 2) return std::sort(begin, end, obj_bb_compare<Item,2> );
+	assert(false);
+}
 
-int rtbuild_size(RTBuilder *b)
+void rtbuild_done(RTBuilder *b)
 {
-	return b->end - b->begin;
+	for(int i=0; i<3; i++)
+	if(b->sorted_begin[i])
+		object_sort( b->sorted_begin[i], b->sorted_end[i], i );
 }
 
+RayObject* rtbuild_get_primitive(RTBuilder *b, int index)
+{
+	return b->sorted_begin[0][index]->obj;
+}
 
 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;
+	rtbuild_init( tmp );
+
+	for(int i=0; i<3; i++)
+		if(b->sorted_begin[i])
+		{
+			tmp->sorted_begin[i] = b->sorted_begin[i] +  b->child_offset[child  ];
+			tmp->sorted_end  [i] = b->sorted_begin[i] +  b->child_offset[child+1];
+		}
+		else
+		{
+			tmp->sorted_begin[i] = 0;
+			tmp->sorted_end  [i] = 0;
+		}
+	
 	return tmp;
 }
 
+void rtbuild_calc_bb(RTBuilder *b)
+{
+	if(b->bb[0] == 1.0e30f)
+	{
+		for(RTBuilder::Object **index = b->sorted_begin[0]; index != b->sorted_end[0]; index++)
+			RE_rayobject_merge_bb( (*index)->obj , 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);
+}
 
 //Left balanced tree
 int rtbuild_mean_split(RTBuilder *b, int nchilds, int axis)
@@ -142,11 +213,13 @@ 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);
@@ -189,120 +262,84 @@ int rtbuild_median_split_largest_axis(RTBuilder *b, int 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];
-};
 
-template<class Obj,int Axis>
-bool obj_bb_compare(const Obj &a, const Obj &b)
-{
-	if(a.bb[Axis] != b.bb[Axis])
-		return a.bb[Axis] < b.bb[Axis];
-	return a.obj < b.obj;
-}
 
-void costobject_sort(CostObject *begin, CostObject *end, int axis)
+struct SweepCost
 {
-	if(axis == 0) return std::sort(begin, end, obj_bb_compare<CostObject,0> );
-	if(axis == 1) return std::sort(begin, end, obj_bb_compare<CostObject,1> );
-	if(axis == 2) return std::sort(begin, end, obj_bb_compare<CostObject,2> );
-	assert(false);
-}
-
-
+	float bb[6];
+	float cost;
+};
 
 /* Object Surface Area Heuristic splitter */
 int rtbuild_heuristic_object_split(RTBuilder *b, int nchilds)
 {
 	int size = rtbuild_size(b);		
 	assert(nchilds == 2);
+	assert(size > 1);
+	int baxis = -1, boffset = 0;
 
-	if(size <= nchilds)
-	{
-		return rtbuild_mean_split_largest_axis(b, nchilds);
-	}
-	else
+	if(size > nchilds)
 	{
 		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" );
+		baxis = -1, boffset = size/2;
 
-		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;
-		}
+		SweepCost *sweep = (SweepCost*)MEM_mallocN( sizeof(SweepCost)*size, "RTBuilder.HeuristicSweep" );
 		
-		for(axis=try_axis[k=0]; k<3; axis=try_axis[++k])
+		for(int axis=0; axis<3; axis++)
 		{
-			float left_cost, right_cost;
-			float other_bb[6];
+			SweepCost sweep_left;
 
-
-
-			costobject_sort(cost, cost+size, axis);
-			for(i=size-1; i>=0; i--)
+			RTBuilder::Object **obj = b->sorted_begin[axis];
+			
+//			float right_cost = 0;
+			for(int 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);
+					VECCOPY(sweep[i].bb, obj[i]->bb);
+					VECCOPY(sweep[i].bb+3, obj[i]->bb+3);
+					sweep[i].cost = obj[i]->cost;
 				}
 				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]);
+					sweep[i].bb[0] = MIN2(obj[i]->bb[0], sweep[i+1].bb[0]);
+					sweep[i].bb[1] = MIN2(obj[i]->bb[1], sweep[i+1].bb[1]);
+					sweep[i].bb[2] = MIN2(obj[i]->bb[2], sweep[i+1].bb[2]);					
+					sweep[i].bb[3] = MAX2(obj[i]->bb[3], sweep[i+1].bb[3]);
+					sweep[i].bb[4] = MAX2(obj[i]->bb[4], sweep[i+1].bb[4]);
+					sweep[i].bb[5] = MAX2(obj[i]->bb[5], sweep[i+1].bb[5]);
+					sweep[i].cost  = obj[i]->cost + sweep[i+1].cost;
 				}
+//				right_cost += obj[i]->cost;
 			}
 			
-			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++)
+			sweep_left.bb[0] = obj[0]->bb[0];
+			sweep_left.bb[1] = obj[0]->bb[1];
+			sweep_left.bb[2] = obj[0]->bb[2];
+			sweep_left.bb[3] = obj[0]->bb[3];
+			sweep_left.bb[4] = obj[0]->bb[4];
+			sweep_left.bb[5] = obj[0]->bb[5];
+			sweep_left.cost  = obj[0]->cost;
+			
+//			right_cost -= obj[0]->cost;	if(right_cost < 0) right_cost = 0;
+
+			for(int 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));
+				left_side = bb_area(sweep_left.bb, sweep_left.bb+3)*(sweep_left.cost+logf(i));
+				right_side= bb_area(sweep[i].bb, sweep[i].bb+3)*(sweep[i].cost+logf(size-i));
+				hcost = left_side+right_side;
+
+				assert(left_side > 0);
+				assert(right_side > 0);
 				
 				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
@@ -311,142 +348,51 @@ int rtbuild_heuristic_object_split(RTBuilder *b, int nchilds)
 					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;
+				DO_MIN( obj[i]->bb,   sweep_left.bb   );
+				DO_MAX( obj[i]->bb+3, sweep_left.bb+3 );
+
+				sweep_left.cost += obj[i]->cost;
+//				right_cost -= obj[i]->cost; if(right_cost < 0) right_cost = 0;
 			}
 			
 			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;		
+		MEM_freeN(sweep);
 	}
-}
-
-/*
- * 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)
+	else if(size == 2)
 	{
-		if(fc < fb)
-			return b;
-		else
-		{
-			if(fc < fa)
-				return c;
-			else
-				return a;
-		}
+		baxis = 0;
+		boffset = 1;
 	}
-	else
+	else if(size == 1)
 	{
-		if(fc < fb)
-		{
-			if(fc < fa)
-				return a;
-			else
-				return c;
-		}
-		else
-			return b;
+		b->child_offset[0] = 0;
+		b->child_offset[1] = 1;
+		return 1;
 	}
-}
-
-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 );
+	b->child_offset[0] = 0;
+	b->child_offset[1] = boffset;
+	b->child_offset[2] = size;
 	
-	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);
+	/* Adjust sorted arrays for childs */
+	for(int i=0; i<boffset; i++) b->sorted_begin[baxis][i]->selected = true;
+	for(int i=boffset; i<size; i++) b->sorted_begin[baxis][i]->selected = false;
+	for(int i=0; i<3; i++)
+		std::stable_partition( b->sorted_begin[i], b->sorted_end[i], selected_node );
 
-	return n;
+	return nchilds;		
 }
 
+/*
+ * Helper code
+ * PARTITION code / used on mean-split
+ * basicly this a std::nth_element (like on C++ STL algorithm)
+ */
+/*
 static void split_leafs(RTBuilder *b, int *nth, int partitions, int split_axis)
 {
 	int i;
@@ -456,23 +402,12 @@ static void split_leafs(RTBuilder *b, int *nth, int partitions, int split_axis)
 	{
 		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++;
-		}
+		if(split_axis == 0)	std::nth_element(b, nth[i],  nth[i+1], nth[partitions], obj_bb_compare<RTBuilder::Object,0>);
+		if(split_axis == 1)	std::nth_element(b, nth[i],  nth[i+1], nth[partitions], obj_bb_compare<RTBuilder::Object,1>);
+		if(split_axis == 2)	std::nth_element(b, nth[i],  nth[i+1], nth[partitions], obj_bb_compare<RTBuilder::Object,2>);
 	}
-	return begin;
 }
+*/
 
 /*
  * Bounding Box utils