# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1195777473 -3600
# Node ID f9fad94cd0cc0d84fe18c2a69d21b6753ef4fe0f
# Parent  3239f749e3940cf59abe2eb9ef38f16f249c52c8
kd-tree: build algorithm tested and fixed
exporting kd-tree to wavefront obj file (visualisation!)

diff -r 3239f749e394 -r f9fad94cd0cc .bzrignore
--- a/.bzrignore	Thu Nov 22 21:46:09 2007 +0100
+++ b/.bzrignore	Fri Nov 23 01:24:33 2007 +0100
@@ -1,1 +1,2 @@
 demos/*.png
+demos/kdtree.obj
diff -r 3239f749e394 -r f9fad94cd0cc src/kdtree.cc
--- a/src/kdtree.cc	Thu Nov 22 21:46:09 2007 +0100
+++ b/src/kdtree.cc	Fri Nov 23 01:24:33 2007 +0100
@@ -22,18 +22,24 @@
 
 void KdNode::subdivide(BBox bbox, int depth)
 {
+	if (depth >= 10 || shapes.size() <= 2)
+	{
+		leaf = true;
+		return;
+	}
+
 	// choose split axis
 	axis = 0;
-	if (bbox.R.y - bbox.L.y > bbox.R.x - bbox.L.x)
+	if (bbox.h() > bbox.w() && bbox.h() > bbox.d())
 		axis = 1;
-	if (bbox.R.z - bbox.L.z > bbox.R.y - bbox.L.y)
+	if (bbox.d() > bbox.w() && bbox.d() > bbox.h())
 		axis = 2;
 
 	// *** find optimal split position
 	SortableShapeList sslist(shapes, axis);
 	sort(sslist.begin(), sslist.end());
 
-	SplitList splitlist = SplitList();
+	SplitList splitlist;
 
 	SortableShapeList::iterator sh;
 	for (sh = sslist.begin(); sh != sslist.end(); sh++)
@@ -42,6 +48,8 @@
 		splitlist.push_back(SplitPos(sh->bbox.R.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;
@@ -51,7 +59,10 @@
 		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.R.cell[axis])
 			{
@@ -69,12 +80,13 @@
 					else
 						spl->lnum++;
 				}
+				sh->setMark();
 			} else
 			// if shape is on left side of split plane
 			if (sh->bbox.R.cell[axis] <= spl->pos)
 			{
+				spl->lnum++;
 				sh->setMark();
-				spl->lnum++;
 			} else
 			// if shape occupies both sides of split plane
 			if (sh->bbox.L.cell[axis] < spl->pos && sh->bbox.R.cell[axis] > spl->pos)
@@ -105,8 +117,8 @@
 		lbb.R.cell[axis] = spl->pos;
 		rbb.L.cell[axis] = spl->pos;
 		float SAL = 2*(lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
-        float SAR = 2*(rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
-        float splitcost = K + SAL/SAV*(K+spl->lnum) + SAR/SAV*(K+spl->rnum);
+		float SAR = 2*(rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
+		float splitcost = K + SAL/SAV*(K+spl->lnum) + SAR/SAV*(K+spl->rnum);
 
 		if (splitcost < cost)
 		{
@@ -119,6 +131,29 @@
 	if (leaf)
 		return;
 
+#if 1
+// export kd-tree as .obj for visualization
+// this would be hard to reconstruct later
+	static ofstream F("kdtree.obj");
+	Vector3 v;
+	static int f=0;
+	v.cell[axis] = splitpos->pos;
+	v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bbox.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.R.cell[(axis+2)%3];
+	F << "v " << v.x << " " << v.y << " " << v.z << endl;
+	v.cell[(axis+1)%3] = bbox.R.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bbox.R.cell[(axis+2)%3];
+	F << "v " << v.x << " " << v.y << " " << v.z << endl;
+	v.cell[(axis+1)%3] = bbox.R.cell[(axis+1)%3];
+	v.cell[(axis+2)%3] = bbox.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;
+#endif
+
 	split = splitpos->pos;
 	float lnum = splitpos->lnum;
 	float rnum = splitpos->rnum;
@@ -184,9 +219,28 @@
 	children[1].subdivide(rbb, depth+1);
 }
 
-void KdTree::build()
+void KdTree::optimize()
 {
 	root = new KdNode();
 	root->shapes = shapes;
 	root->subdivide(bbox, 0);
+	built = true;
 }
+
+void KdTree::save(ostream &str, KdNode *node)
+{
+	if (!built)
+		return;
+	if (node == NULL)
+		node = root;
+	if (node->isLeaf())
+		str << "(leaf: " << node->shapes.size() << " shapes)";
+	else
+	{
+		str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L=";
+		save(str, node->getLeftChild());
+		str << "; R=";
+		save(str, node->getRightChild());
+		str << ";)";
+	}
+}
diff -r 3239f749e394 -r f9fad94cd0cc src/kdtree.h
--- a/src/kdtree.h	Thu Nov 22 21:46:09 2007 +0100
+++ b/src/kdtree.h	Fri Nov 23 01:24:33 2007 +0100
@@ -2,6 +2,8 @@
 #define KDTREE_H
 
 #include <vector>
+#include <iostream>
+#include <fstream>
 
 #include "scene.h"
 
@@ -44,9 +46,12 @@
 public:
 	float pos;
 	int lnum, rnum;
-	SplitPos(float &aPos): pos(aPos) {};
+	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>
@@ -56,12 +61,13 @@
 class Container
 {
 protected:
-	ShapeList shapes;
 	BBox bbox;
 public:
-	Container(): shapes(), bbox() {};
-	void addShape(Shape* aShape);
+	ShapeList shapes;
+	Container(): bbox(), shapes() {};
+	virtual void addShape(Shape* aShape);
 	//void addShapeNoExtend(Shape* aShape) { shapes.push_back(aShape); };
+	virtual void optimize() {};
 };
 
 class KdNode
@@ -95,9 +101,12 @@
 class KdTree: public Container
 {
 	KdNode *root;
+	bool built;
 public:
-	KdTree() {};
-	void build();
+	KdTree() : Container(), built(false) {};
+	void addShape(Shape* aShape) { Container::addShape(aShape); built = false; };
+	void optimize(); // build kd-tree
+	void save(ostream &str, KdNode *node = NULL);
 };
 
 #endif
diff -r 3239f749e394 -r f9fad94cd0cc src/raytracer.cc
--- a/src/raytracer.cc	Thu Nov 22 21:46:09 2007 +0100
+++ b/src/raytracer.cc	Fri Nov 23 01:24:33 2007 +0100
@@ -57,7 +57,7 @@
 {
 	Shape *nearest_shape = NULL;
 	ShapeList::iterator shape;
-	for (shape = shapes.begin(); shape != shapes.end(); shape++)
+	for (shape = top->shapes.begin(); shape != top->shapes.end(); shape++)
 		if (*shape != origin_shape && (*shape)->intersect(ray, nearest_distance))
 			nearest_shape = *shape;
 	return nearest_shape;
@@ -221,6 +221,12 @@
 	float vy;
 	RenderrowData *d;
 
+	printf("* building kd-tree\n");
+	//cout << endl;
+	static_cast<KdTree*>(top)->optimize();
+	//static_cast<KdTree*>(top)->save(cout);
+	//cout << endl;
+
 	//srand(time(NULL));
 
 	// eye - viewing point
@@ -283,11 +289,6 @@
 	return data;
 }
 
-void Raytracer::addshape(Shape *shape)
-{
-	shapes.push_back(shape);
-}
-
 void Raytracer::addlight(Light *light)
 {
 	lights.push_back(light);
diff -r 3239f749e394 -r f9fad94cd0cc src/raytracer.h
--- a/src/raytracer.h	Thu Nov 22 21:46:09 2007 +0100
+++ b/src/raytracer.h	Fri Nov 23 01:24:33 2007 +0100
@@ -27,7 +27,7 @@
 
 class Raytracer
 {
-	ShapeList shapes;
+	Container *top;
 	vector<Light*> lights;
 	Colour bg_colour;
 	int ao_samples;
@@ -37,11 +37,11 @@
 	inline Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
 		float &nearest_distance);
 public:
-	Raytracer(): shapes(), lights(), bg_colour(0.0, 0.0, 0.0), ao_samples(0) {};
+	Raytracer(): lights(), bg_colour(0.0, 0.0, 0.0), ao_samples(0) { top = new KdTree(); };
 	~Raytracer() {};
 	float *render(int w, int h);
 	Colour raytrace(Ray &ray, int depth, Shape *origin_shape);
-	void addshape(Shape *shape);
+	void addshape(Shape *shape) { top->addShape(shape); };
 	void addlight(Light *light);
 
 	void ambientocclusion(int samples, float distance, float angle);