# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1208620827 -7200
# Node ID 4fedf72909297cafa31ac3d8dd13afdb799b5449
# Parent  4b84e90325c5d3a665f6e1b9112017474cc3301a
simplify kd-tree building, it's also much faster now

diff -r 4b84e90325c5 -r 4fedf7290929 ccdemos/realtime_dragon.cc
--- a/ccdemos/realtime_dragon.cc	Tue Apr 15 17:12:50 2008 +0200
+++ b/ccdemos/realtime_dragon.cc	Sat Apr 19 18:00:27 2008 +0200
@@ -1,5 +1,5 @@
 #include "raytracer.h"
-#include "octree.h"
+#include "kdtree.h"
 
 #include "common_sdl.h"
 #include "common_ply.h"
@@ -7,7 +7,7 @@
 int main(int argc, char **argv)
 {
 	Raytracer rt;
-	Octree top;
+	KdTree top;
 	Camera cam;
 
 	rt.setMaxDepth(0);
@@ -30,5 +30,8 @@
 
 	top.optimize();
 
+	if (argc == 2 && strcmp(argv[1], "-buildonly") == 0)
+		return 0;
+
 	loop_sdl(rt, cam);
 }
diff -r 4b84e90325c5 -r 4fedf7290929 src/kdtree.cc
--- a/src/kdtree.cc	Tue Apr 15 17:12:50 2008 +0200
+++ b/src/kdtree.cc	Sat Apr 19 18:00:27 2008 +0200
@@ -30,51 +30,23 @@
 #include "common.h"
 #include "kdtree.h"
 
-class SortableShape
+class ShapeBound
 {
 public:
 	Shape *shape;
-	BBox bbox;
-	short axis;
-	short mark;
-
-	SortableShape(Shape *aShape, short &aAxis): shape(aShape), axis(aAxis), mark(0)
-		{ bbox = shape->get_bbox(); };
-	friend bool operator<(const SortableShape& a, const SortableShape& b)
-		{ return a.bbox.L.cell[a.axis] < b.bbox.L.cell[b.axis]; };
-	void setAxis(short aAxis) { axis = aAxis; };
-	void setMark() { mark = 1; };
-	short hasMark() { return mark; };
-};
-
-class SortableShapeList: public vector<SortableShape>
-{
-public:
-	SortableShapeList(ShapeList *shapes, short axis)
+	Float pos;
+	bool end;
+	ShapeBound(Shape *ashape, const Float apos, const bool aend):
+		shape(ashape), pos(apos), end(aend) {};
+	friend bool operator<(const ShapeBound& a, const ShapeBound& b)
 	{
-		ShapeList::iterator shape;
-		for (shape = shapes->begin(); shape != shapes->end(); shape++)
-			push_back(SortableShape(*shape, axis));
+		if (a.pos == b.pos)
+			return a.shape < b.shape;
+		else
+			return a.pos < b.pos;
 	};
 };
 
-class SplitPos
-{
-public:
-	Float pos;
-	int lnum, rnum;
-	SplitPos(): pos(0.0), lnum(0), rnum(0) {};
-	SplitPos(Float &aPos): pos(aPos), lnum(0), rnum(0) {};
-	friend bool operator<(const SplitPos& a, const SplitPos& b)
-		{ return a.pos < b.pos; };
-	friend bool operator==(const SplitPos& a, const SplitPos& b)
-		{ return a.pos == b.pos; };
-};
-
-class SplitList: public vector<SplitPos>
-{
-};
-
 // stack element for kd-tree traversal
 class StackElem
 {
@@ -96,7 +68,8 @@
 		delete[] children;
 }
 
-void KdNode::subdivide(BBox bbox, int maxdepth)
+// kd-tree recursive build algorithm, inspired by PBRT (www.pbrt.org)
+void KdNode::subdivide(BBox bounds, int maxdepth)
 {
 	if (maxdepth <= 0 || shapes->size() <= 2)
 	{
@@ -106,104 +79,59 @@
 
 	// choose split axis
 	axis = 0;
-	if (bbox.h() > bbox.w() && bbox.h() > bbox.d())
+	if (bounds.h() > bounds.w() && bounds.h() > bounds.d())
 		axis = 1;
-	if (bbox.d() > bbox.w() && bbox.d() > bbox.h())
+	if (bounds.d() > bounds.w() && bounds.d() > bounds.h())
 		axis = 2;
 
-	// *** find optimal split position
-	SortableShapeList sslist(shapes, axis);
-	sort(sslist.begin(), sslist.end());
-
-	SplitList splitlist;
-
-	SortableShapeList::iterator sh;
-	for (sh = sslist.begin(); sh != sslist.end(); sh++)
-	{
-		splitlist.push_back(SplitPos(sh->bbox.L.cell[axis]));
-		splitlist.push_back(SplitPos(sh->bbox.H.cell[axis]));
-	}
-	sort(splitlist.begin(), splitlist.end());
-	// remove duplicate splits
-	splitlist.resize(unique(splitlist.begin(), splitlist.end()) - splitlist.begin());
-
-	// find all posible splits
-	SplitList::iterator spl;
-	int rest;
-	for (spl = splitlist.begin(); spl != splitlist.end(); spl++)
+	// create sorted list of shape bounds (= find all posible splits)
+	vector<ShapeBound> edges[3];
+	ShapeList::iterator shape;
+	for (shape = shapes->begin(); shape != shapes->end(); shape++)
 	{
-		for (sh = sslist.begin(), rest = sslist.size(); sh != sslist.end(); sh++, rest--)
-		{
-			if (sh->hasMark())
-			{
-				spl->lnum++;
-				continue;
-			}
-			// if shape is completely contained in split plane
-			if (spl->pos == sh->bbox.L.cell[axis] == sh->bbox.H.cell[axis])
-			{
-				if (spl->pos - bbox.L.cell[axis] < bbox.H.cell[axis] - spl->pos)
-				{
-					// left subcell is smaller -> if not empty, put shape here
-					if (spl->lnum)
-						spl->lnum++;
-					else
-						spl->rnum++;
-				} else {
-					// right subcell is smaller
-					if (spl->rnum)
-						spl->rnum++;
-					else
-						spl->lnum++;
-				}
-				sh->setMark();
-			} else
-			// if shape is on left side of split plane
-			if (sh->bbox.H.cell[axis] <= spl->pos)
-			{
-				spl->lnum++;
-				sh->setMark();
-			} else
-			// if shape occupies both sides of split plane
-			if (sh->bbox.L.cell[axis] < spl->pos && sh->bbox.H.cell[axis] > spl->pos)
-			{
-				spl->lnum++;
-				spl->rnum++;
-			} else
-			// if shape is on right side of split plane
-			if (sh->bbox.L.cell[axis] >= spl->pos)
-			{
-				spl->rnum += rest;
-				break;
-			}
-		}
+		BBox shapebounds = (*shape)->get_bbox();
+	//	for (int ax = 0; ax < 3; ax++)
+	//	{
+			edges[axis].push_back(ShapeBound(*shape, shapebounds.L[axis], 0));
+			edges[axis].push_back(ShapeBound(*shape, shapebounds.H[axis], 1));
+	//	}
 	}
+	sort(edges[axis].begin(), edges[axis].end());
 
 	// choose best split pos
 	const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
-	Float SAV = (bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node
+	Float SAV = (bounds.w()*bounds.h() +  // surface area of node
+		bounds.w()*bounds.d() + bounds.h()*bounds.d());
 	Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost
-	bool leaf = true;
-	BBox lbb = bbox;
-	BBox rbb = bbox;
-	for (spl = splitlist.begin(); spl != splitlist.end(); spl++)
+	BBox lbb = bounds;
+	BBox rbb = bounds;
+
+	vector<ShapeBound>::iterator edge, splitedge = edges[axis].end();
+	int lnum = 0, rnum = shapes->size();
+	for (edge = edges[axis].begin(); edge != edges[axis].end(); edge++)
 	{
+		if (edge->end)
+			rnum--;
+
 		// calculate SAH cost of this split
-		lbb.H.cell[axis] = spl->pos;
-		rbb.L.cell[axis] = spl->pos;
+		lbb.H.cell[axis] = edge->pos;
+		rbb.L.cell[axis] = edge->pos;
 		Float SAL = (lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
 		Float SAR = (rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
-		Float splitcost = K*SAV + SAL*(K+spl->lnum) + SAR*(K+spl->rnum);
+		Float splitcost = K*SAV + SAL*(K + lnum) + SAR*(K + rnum);
 
 		if (splitcost < cost)
 		{
-			leaf = false;
+			splitedge = edge;
 			cost = splitcost;
-			split = spl->pos;
+			split = edge->pos;
 		}
+
+		if (!edge->end)
+			lnum++;
 	}
 
-	if (leaf)
+	if (splitedge == edges[axis].end())
 	{
 		setLeaf();
 		return;
@@ -216,17 +144,17 @@
 	Vector3 v;
 	static int f=0;
 	v.cell[axis] = split;
-	v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3];
-	v.cell[(axis+2)%3] = bbox.L.cell[(axis+2)%3];
+	v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bounds.L.cell[(axis+2)%3];
 	F << "v " << v.x << " " << v.y << " " << v.z << endl;
-	v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3];
-	v.cell[(axis+2)%3] = bbox.H.cell[(axis+2)%3];
+	v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bounds.H.cell[(axis+2)%3];
 	F << "v " << v.x << " " << v.y << " " << v.z << endl;
-	v.cell[(axis+1)%3] = bbox.H.cell[(axis+1)%3];
-	v.cell[(axis+2)%3] = bbox.H.cell[(axis+2)%3];
+	v.cell[(axis+1)%3] = bounds.H.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bounds.H.cell[(axis+2)%3];
 	F << "v " << v.x << " " << v.y << " " << v.z << endl;
-	v.cell[(axis+1)%3] = bbox.H.cell[(axis+1)%3];
-	v.cell[(axis+2)%3] = bbox.L.cell[(axis+2)%3];
+	v.cell[(axis+1)%3] = bounds.H.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bounds.L.cell[(axis+2)%3];
 	F << "v " << v.x << " " << v.y << " " << v.z << endl;
 	F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl;
 	f += 4;
@@ -235,57 +163,12 @@
 	// split this node
 	delete shapes;
 	children = new KdNode[2];
-	int state = 0;
-	for (sh = sslist.begin(); sh != sslist.end(); sh++)
-	{
-		// if shape is on left side of split plane
-		if (state == 1)
-		{ // only right
-			children[1].addShape(sh->shape);
-			continue;
-		}
-		if (state == 0)
-		{
-			if (sh->bbox.H.cell[axis] < split)
-			{ // left
-				children[0].addShape(sh->shape);
-			} else
-			if (sh->bbox.H.cell[axis] > split)
-			{
-				if (sh->bbox.L.cell[axis] < split)
-				{ // both
-					children[0].addShape(sh->shape);
-					children[1].addShape(sh->shape);
-				} else
-				{ // right
-					children[1].addShape(sh->shape);
-					state = 1;
-				}
-			} else
-			{ // H == split
-				if (sh->bbox.L.cell[axis] < split)
-				{ // left
-					children[0].addShape(sh->shape);
-				} else
-				{ // contained
-					if (split - bbox.L.cell[axis] < bbox.H.cell[axis] - split)
-					{
-						// left subcell is smaller -> if not empty, put shape here
-						if (children[0].shapes->size())
-							children[0].addShape(sh->shape);
-						else
-							children[1].addShape(sh->shape);
-					} else {
-						// right subcell is smaller
-						if (children[1].shapes->size())
-							children[1].addShape(sh->shape);
-						else
-							children[0].addShape(sh->shape);
-					}
-				}
-			}
-		}
-	}
+	for (edge = edges[axis].begin(); edge != splitedge; edge++)
+		if (!edge->end)
+			children[0].addShape(edge->shape);
+	for (edge = splitedge; edge < edges[axis].end(); edge++)
+		if (edge->end)
+			children[1].addShape(edge->shape);
 
 	lbb.H.cell[axis] = split;
 	rbb.L.cell[axis] = split;