# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1195750414 -3600
# Node ID bf17f9f84c916ffb4fac1959ef6924573b3e5379
# Parent  d8d596d26f25273a70dbc17b8279a9b90c84ef54
kd-tree: build algorithm - searching for all posible splits

diff -r d8d596d26f25 -r bf17f9f84c91 DEVNOTES
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DEVNOTES	Thu Nov 22 17:53:34 2007 +0100
@@ -0,0 +1,58 @@
+KdTree Algorithm
+----------------
+
+def build_kdtree(root, shapes):
+	root = new KdNode
+	root.shapes = shapes
+	subdivide(root, bounding_box(shapes), 0)
+
+def subdivide(node, bbox, depth):
+	# choose split axis
+	axis = bbox.largest_extend()
+
+	# find split position
+	splitpos = find_split_position(axis)
+
+	leftnode, rightnode = new KdNode[2]
+	for each shape:
+		if (node->intersectleftnode()) leftnode->addprimitive( primitive )
+		if (node->intersectrightnode()) rightnode->addprimitive( primitive )
+	subdivide(leftnode, bbox.left_split(splitpos), depth+1)
+	subdivide(rightnode, bbox.right_split(splitpos), depth+1)
+
+def find_split_position(axis):
+	mshapes = new ShapeListWithMarks(shapes)
+
+	# sort new shape list according to left boundaries
+	mshapes.sort(axis)
+
+	splitlist = new SplitList()
+	for each mshape from mshapes:
+		splitlist.add_extremes(mshape, axis)
+	splitlist.sort()
+	for each split from splitlist:
+		for each mshape from mshapes:
+			if mshape.marked:
+				continue
+
+			# if shape is on left side of split plane
+			if mshape.r_boundary <= split.pos:
+				split.lnum++
+				mshape.mark()
+
+			# if shape is completely contained in split plane
+			if mshape.l_boundary == mshape.r_boundary == split.pos:
+				if split.num_smaller_subcell == 0:
+					split.num_bigger_subcell++
+				else:
+					split.num_smaller_subcell++
+
+			# if shape is on right side of split plane
+			if mshape.l_boundary >= split.pos:
+				split.r_num += mshapes.count - mshape.number
+				break
+
+			# if shape occupies both sides of split plane
+			if mshape.l_boundary < split.pos and mshape.r_boundary > split.pos:
+				split.l_num++
+				split.r_num++
diff -r d8d596d26f25 -r bf17f9f84c91 Makefile
--- a/Makefile	Sun Nov 18 11:20:56 2007 +0100
+++ b/Makefile	Thu Nov 22 17:53:34 2007 +0100
@@ -46,11 +46,12 @@
 matrix.o: src/matrix.cc src/matrix.h src/vector.h
 noise.o: src/noise.cc src/noise.h
 scene.o: src/scene.cc src/scene.h src/vector.h src/noise.h
+kdtree.o: src/kdtree.cc src/kdtree.h
 raytracer.o: src/raytracer.cc src/raytracer.h src/scene.h src/vector.h src/noise.h
 
 # python module
 raytracermodule.o: src/raytracermodule.cc src/raytracer.h src/scene.h src/vector.h
-$(MODULENAME): raytracermodule.o raytracer.o scene.o noise.o
+$(MODULENAME): raytracermodule.o raytracer.o scene.o noise.o kdtree.o
 	$(CXX) $^ -shared -o $@ $(LDFLAGS)
 
 # library tests
diff -r d8d596d26f25 -r bf17f9f84c91 TODO
--- a/TODO	Sun Nov 18 11:20:56 2007 +0100
+++ b/TODO	Thu Nov 22 17:53:34 2007 +0100
@@ -1,5 +1,8 @@
  * kd-tree
+ * transforms, camera
 
+New Classes?
+============
 
 container.h  -- Container
 kdtree.h  -- KdTree
@@ -10,17 +13,19 @@
 vector.h  -- Vector3
 reader.h  -- Reader, WavefrontReader
 
-KdTree monkey
+KdTree top
 wf = new WavefrontReader()
-wf.setContainer(monkey)
+wf.setContainer(top)
+wf.setTransform(monkey_pos_matrix)
 wf.read("monkey.obj")
+// more transform&reads
 destroy wf
 
-monkey.optimize()  -- i.e. build tree
+top.optimize()  -- i.e. build tree
 
 Scene scene -- container with shapes, a camera and lights
 scene = new Scene()
-scene.setTop(monkey) -- top object in hierarchy
+scene.setTop(top) -- top object in hierarchy
 scene.setCamera(new Camera(pos, dir, angle))
 scene.addLight(new PointLight(pos, color))
 rt.setScene(scene)
diff -r d8d596d26f25 -r bf17f9f84c91 src/kdtree.cc
--- a/src/kdtree.cc	Sun Nov 18 11:20:56 2007 +0100
+++ b/src/kdtree.cc	Thu Nov 22 17:53:34 2007 +0100
@@ -1,36 +1,129 @@
-#include "kdtree.cc"
+#include <algorithm>
+
+#include "kdtree.h"
 
-void KdTree::KdTree(ShapeList &shapelist):
+void Container::addShape(Shape* aShape)
+{
+	shapes.push_back(aShape);
+	if (shapes.size() == 0) {
+		/* initialize bounding box */
+		bbox = aShape->get_bbox();
+	} else {
+		/* adjust bounding box */
+		BBox shapebb = aShape->get_bbox();
+		if (shapebb.L.x < bbox.L.x)  bbox.L.x = shapebb.L.x;
+		if (shapebb.L.y < bbox.L.y)  bbox.L.y = shapebb.L.y;
+		if (shapebb.L.z < bbox.L.z)  bbox.L.z = shapebb.L.z;
+		if (shapebb.R.x > bbox.R.x)  bbox.R.x = shapebb.R.x;
+		if (shapebb.R.y > bbox.R.y)  bbox.R.y = shapebb.R.y;
+		if (shapebb.R.z > bbox.R.z)  bbox.R.z = shapebb.R.z;
+	}
+};
+
+void KdNode::subdivide(BBox bbox, int depth, int count)
 {
-	root = new KdNode();
-	shapes = new vector<Shape*>();
+	int axis;
+	float pos;
+
+	//if (stopcriterionmet()) return
+
+	// choose split axis
+	axis = 0;
+	if (bbox.R.y - bbox.L.y > bbox.R.x - bbox.L.x)
+		axis = 1;
+	if (bbox.R.z - bbox.L.z > bbox.R.y - bbox.L.y)
+		axis = 2;
+
+	// *** find optimal split position
+	SortableShapeList sslist(shapes, axis);
+	sort(sslist.begin(), sslist.end());
+
+	SplitList 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.R.cell[axis]));
+	}
+	sort(splitlist.begin(), splitlist.end());
+
+	// find all posible splits
+	SplitList::iterator split;
+	int rest;
+	for (split = splitlist.begin(); split != splitlist.end(); split++)
+	{
+		for (sh = sslist.begin(), rest = sslist.size(); sh != sslist.end(); sh++, rest--)
+		{
+			if (sh->hasMark())
+				continue;
+			// if shape is on left side of split plane
+			if (sh->bbox.R.cell[axis] <= split->pos)
+			{
+				sh->setMark();
+				split->lnum++;
+			}
+			// if shape is completely contained in split plane
+			if (sh->bbox.L.cell[axis] == sh->bbox.R.cell[axis] == split->pos)
+			{
+				if (split->pos - bbox.L.cell[axis] < bbox.R.cell[axis] - split->pos)
+				{
+					// left subcell is smaller -> if not empty, put shape here
+					if (split->lnum)
+						split->lnum++;
+					else
+						split->rnum++;
+				} else {
+					// right subcell is smaller
+					if (split->rnum)
+						split->rnum++;
+					else
+						split->lnum++;
+				}
+			}
+			// if shape is on right side of split plane
+			if (sh->bbox.L.cell[axis] >= split->pos)
+			{
+				split->rnum += rest;
+				break;
+			}
+			// if shape occupies both sides of split plane
+			if (sh->bbox.L.cell[axis] < split->pos && sh->bbox.R.cell[axis] > split->pos)
+			{
+				split->lnum++;
+				split->rnum++;
+			}
+		}
+	}
+
+	// choose best split pos
+	// ... //
+
+	KdNode *children = new KdNode[2];
+
 	ShapeList::iterator shape;
-	for (shape = shapelist.begin(); shape != shapes.end(); shape++)
-		addShape(*shape);
+	for (shape = shapes.begin(); shape != shapes.end(); shape++)
+	{
+		if (((Triangle*)*shape)->A.cell[axis-1] < pos
+		|| ((Triangle*)*shape)->B.cell[axis-1] < pos
+		|| ((Triangle*)*shape)->C.cell[axis-1] < pos)
+			children[0].addShape(*shape);
+		if (((Triangle*)*shape)->A.cell[axis-1] >= pos
+		|| ((Triangle*)*shape)->B.cell[axis-1] >= pos
+		|| ((Triangle*)*shape)->C.cell[axis-1] >= pos)
+			children[1].addShape(*shape);
+	}
 
-	rebuild();
+	bbox.R[axis-1] = pos;
+	children[0].subdivide(bbox, depth+1, children[0].shapes.size());
+
+	bbox.L[axis-1] = pos;
+	children[1].subdivide(bbox, depth+1, children[1].shapes.size());
 }
 
-void KdTree::Subdivide(KdNode* node, AABB& bbox, int depth, int count)
+void KdTree::build()
 {
-	
-	/*if (stopcriterionmet()) return
-	splitpos = findoptimalsplitposition()
-	leftnode = new Node()
-	rightnode = new Node()
-	for (all primitives in node)
-	{
-		if (node->intersectleftnode()) leftnode->addprimitive( primitive )
-		if (node->intersectrightnode()) rightnode->addprimitive( primitive )
-	}
-	buildkdtree( leftnode )
-		buildkdtree( rightnode )*/
+	root = new KdNode();
+	root->shapes = shapes;
+	root->subdivide(bbox, 0, shapes.size());
 }
-
-void KdTree::rebuild()
-{
-	int count = shapes->size();
-	AABB bbox = shapes->extends();
-
-	Subdivide(root, bbox, 0, count);
-}
\ No newline at end of file
diff -r d8d596d26f25 -r bf17f9f84c91 src/kdtree.h
--- a/src/kdtree.h	Sun Nov 18 11:20:56 2007 +0100
+++ b/src/kdtree.h	Thu Nov 22 17:53:34 2007 +0100
@@ -5,47 +5,100 @@
 
 #include "scene.h"
 
-class SpaceDivider
+using namespace std;
+
+class ShapeList: public vector<Shape*>
+{
+};
+
+class SortableShape
 {
-	ShapeList *shapes;
 public:
-	SpaceDivider(ShapeList &shapelist): shapes(shapelist) {};
+	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 KdNode:
+class SortableShapeList: public vector<SortableShape>
+{
+public:
+	SortableShapeList(ShapeList &shapes, short axis)
+	{
+		ShapeList::iterator shape;
+		for (shape = shapes.begin(); shape != shapes.end(); shape++)
+			push_back(SortableShape(*shape, axis));
+	};
+};
+
+class SplitPos
+{
+public:
+	float pos;
+	int lnum, rnum;
+	SplitPos(float &aPos): pos(aPos) {};
+	friend bool operator<(const SplitPos& a, const SplitPos& b)
+		{ return a.pos < b.pos; };
+};
+
+class SplitList: public vector<SplitPos>
+{
+};
+
+class Container
+{
+protected:
+	ShapeList shapes;
+	BBox bbox;
+public:
+	Container(): shapes(), bbox() {};
+	void addShape(Shape* aShape);
+	//void addShapeNoExtend(Shape* aShape) { shapes.push_back(aShape); };
+};
+
+class KdNode
 {
 	float split;
 	bool leaf; /* is this node a leaf? */
-	char axis; /* 0,1,2 => x,y,z */
-	KdNode *leftchild, *rightchild;
+	char axis; /* 1,2,3 => x,y,z */
+	KdNode *leftchild;
 public:
-	vector<Shape*> shapes;
+	ShapeList shapes;
 
 	KdNode() : leaf(true), axis(0), shapes() {};
 
-	setAxis(char aAxis) { axis = aAxis; };
+	void setAxis(char aAxis) { axis = aAxis; };
 	char getAxis() { return axis; };
 
-	setSplit(float aSplit) { split = aSplit; };
+	void setSplit(float aSplit) { split = aSplit; };
 	float getSplit() { return split; };
 
-	setLeaf(bool aLeaf) { leaf = aLeaf; };
+	void setLeaf(bool aLeaf) { leaf = aLeaf; };
 	bool isLeaf() { return leaf; };
 
-	setLeftChild(KdNode *aLeft) { leftchild = aLeft; };
+	void setLeftChild(KdNode *aLeft) { leftchild = aLeft; };
 	KdNode *getLeftChild() { return leftchild; };
-	setRightChild(KdNode *aRight) { rightchild = aRight; };
-	KdNode *getRightChild() { return rightchild; };
+	KdNode *getRightChild() { return leftchild+1; };
 
-	addShape(Shape* aShape) { shapes.push_back(aShape); };
+	void addShape(Shape* aShape) { shapes.push_back(aShape); };
+
+	void subdivide(BBox bbox, int depth, int count);
 };
 
-class KdTree: public SpaceDivider
+class KdTree: public Container
 {
-	KdNote *root;
+	KdNode *root;
 public:
-	KdTree(ShapeList &shapelist);
-	rebuild();
+	KdTree() {};
+	void build();
 };
 
 #endif
diff -r d8d596d26f25 -r bf17f9f84c91 src/raytracer.h
--- a/src/raytracer.h	Sun Nov 18 11:20:56 2007 +0100
+++ b/src/raytracer.h	Thu Nov 22 17:53:34 2007 +0100
@@ -10,20 +10,11 @@
 
 #include <vector>
 
+#include "kdtree.h"
 #include "scene.h"
 
 using namespace std;
 
-/* axis-aligned bounding box */
-class AABB
-{
-};
-
-class ShapeList: public vector<Shape*>
-{
-	AABB extends() { return AABB(); };
-};
-
 class Raytracer;
 struct RenderrowData {
 	Raytracer *rt;
diff -r d8d596d26f25 -r bf17f9f84c91 src/scene.cc
--- a/src/scene.cc	Sun Nov 18 11:20:56 2007 +0100
+++ b/src/scene.cc	Thu Nov 22 17:53:34 2007 +0100
@@ -165,3 +165,23 @@
 	dist = t;
 	return true;
 }
+
+BBox Triangle::get_bbox()
+{
+	BBox bbox = BBox();
+	bbox.L = A;
+	if (B.x < bbox.L.x)  bbox.L.x = B.x;
+	if (C.x < bbox.L.x)  bbox.L.x = C.x;
+	if (B.y < bbox.L.y)  bbox.L.y = B.y;
+	if (C.y < bbox.L.y)  bbox.L.y = C.y;
+	if (B.z < bbox.L.z)  bbox.L.z = B.z;
+	if (C.z < bbox.L.z)  bbox.L.z = C.z;
+	bbox.R = A;
+	if (B.x > bbox.R.x)  bbox.R.x = B.x;
+	if (C.x > bbox.R.x)  bbox.R.x = C.x;
+	if (B.y > bbox.R.y)  bbox.R.y = B.y;
+	if (C.y > bbox.R.y)  bbox.R.y = C.y;
+	if (B.z > bbox.R.z)  bbox.R.z = B.z;
+	if (C.z > bbox.R.z)  bbox.R.z = C.z;
+	return bbox;
+};
diff -r d8d596d26f25 -r bf17f9f84c91 src/scene.h
--- a/src/scene.h	Sun Nov 18 11:20:56 2007 +0100
+++ b/src/scene.h	Thu Nov 22 17:53:34 2007 +0100
@@ -16,6 +16,15 @@
 
 using namespace std;
 
+/* axis-aligned bounding box */
+class BBox
+{
+public:
+	Vector3 L;
+	Vector3 R;
+	BBox(): L(), R() {};
+};
+
 class Ray
 {
 public:
@@ -84,6 +93,8 @@
 
 	// normal at point P
 	virtual Vector3 normal(Vector3 &P) = 0;
+
+	virtual BBox get_bbox();
 };
 
 class Sphere: public Shape
@@ -130,6 +141,7 @@
 	bool intersect(const Ray &ray, float &dist);
 	bool intersect_all(const Ray &ray, float dist, vector<float> &allts) {return false;};
 	Vector3 normal(Vector3 &) { return N; };
+	BBox get_bbox();
 };
 
 #endif