# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1197485959 -3600
# Node ID fb170fccb19f10f10712611e92a02ee436daa2a3
# Parent  28f6e8b9d5d1f68172caa9b86279e5cdfbccfe83
new space partitioning structure: octree
realtime_bunny updated to use octree
plus other files updated to be container type independent (only user programs are supposed to include and use special containers)

diff -r 28f6e8b9d5d1 -r fb170fccb19f ccdemos/realtime_bunny.cc
--- a/ccdemos/realtime_bunny.cc	Sun Dec 09 23:10:58 2007 +0100
+++ b/ccdemos/realtime_bunny.cc	Wed Dec 12 19:59:19 2007 +0100
@@ -1,6 +1,7 @@
 #include <SDL.h>
 
 #include "raytracer.h"
+#include "octree.h"
 #include <iostream>
 #include <fstream>
 #include <iomanip>
@@ -144,10 +145,10 @@
 	/* initialize raytracer and prepare scene */
 	render_buffer = (Float *) malloc(w*h*3*sizeof(Float));
 
-	rt.setThreads(1);
+	rt.setThreads(2);
 	rt.setMaxDepth(3);
 
-	KdTree top;
+	Octree top;
 	rt.setTop(&top);
 
 	Light light1(Vector3(-5.0, 2.0, 8.0), Colour(0.9, 0.3, 0.6));
@@ -165,7 +166,7 @@
 	cam.setEye(Vector3(0,0,10));
 
 	/* build kd-tree */
-	top.setMaxDepth(100);
+	top.setMaxDepth(3);
 	top.optimize();
 
 	/* loop... */
diff -r 28f6e8b9d5d1 -r fb170fccb19f ccdemos/realtime_dragon.cc
--- a/ccdemos/realtime_dragon.cc	Sun Dec 09 23:10:58 2007 +0100
+++ b/ccdemos/realtime_dragon.cc	Wed Dec 12 19:59:19 2007 +0100
@@ -1,6 +1,7 @@
 #include <SDL.h>
 
 #include "raytracer.h"
+#include "kdtree.h"
 #include <iostream>
 #include <fstream>
 
@@ -142,7 +143,7 @@
 	rt.setThreads(1);
 	rt.setMaxDepth(3);
 
-	KdTree top;
+	Kdtree top;
 	rt.setTop(&top);
 
 	Light light1(Vector3(-5.0, 2.0, 8.0), Colour(0.9, 0.3, 0.6));
@@ -160,7 +161,7 @@
 	cam.setEye(Vector3(0,0,10));
 
 	/* build kd-tree */
-	top.setMaxDepth(100);
+	top.setMaxDepth(30);
 	top.optimize();
 
 	/* loop... */
diff -r 28f6e8b9d5d1 -r fb170fccb19f include/kdtree.h
--- a/include/kdtree.h	Sun Dec 09 23:10:58 2007 +0100
+++ b/include/kdtree.h	Wed Dec 12 19:59:19 2007 +0100
@@ -44,7 +44,7 @@
 
 	void addShape(Shape* aShape) { shapes->push_back(aShape); };
 
-	void subdivide(BBox bbox, int depth);
+	void subdivide(BBox bbox, int maxdepth);
 };
 
 class KdTree: public Container
diff -r 28f6e8b9d5d1 -r fb170fccb19f include/octree.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/octree.h	Wed Dec 12 19:59:19 2007 +0100
@@ -0,0 +1,53 @@
+/*
+ * Pyrit Ray Tracer
+ * file: octree.h
+ *
+ * Radek Brich, 2006-2007
+ */
+
+#ifndef OCTREE_H
+#define OCTREE_H
+
+#include "container.h"
+#include "vector.h"
+#include "scene.h"
+
+using namespace std;
+
+class OctreeNode
+{
+	OctreeNode *children; // pointer to first of eight children
+public:
+	ShapeList *shapes;    // pointer to shape array, if this is leaf
+
+	OctreeNode() : children(NULL) { shapes = new ShapeList(); };
+	~OctreeNode();
+
+	bool isLeaf() { return shapes != NULL; };
+
+	OctreeNode *getChild(const int num) { return children+num; };
+
+	void addShape(Shape* aShape) { shapes->push_back(aShape); };
+
+	void subdivide(BBox bbox, int maxdepth);
+};
+
+class Octree: public Container
+{
+	OctreeNode *root;
+	bool built;
+	int max_depth;
+public:
+	Octree() : Container(), root(NULL), built(false), max_depth(4) {};
+	~Octree() { if (root) delete root; };
+	void addShape(Shape* aShape) { Container::addShape(aShape); built = false; };
+	Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
+		Float &nearest_distance);
+	void optimize() { build(); };
+	void build();
+	void save(ostream &str, OctreeNode *node = NULL) {};
+	void load(istream &str, OctreeNode *node = NULL) {};
+	void setMaxDepth(int md) { max_depth = md; };
+};
+
+#endif
diff -r 28f6e8b9d5d1 -r fb170fccb19f include/raytracer.h
--- a/include/raytracer.h	Sun Dec 09 23:10:58 2007 +0100
+++ b/include/raytracer.h	Wed Dec 12 19:59:19 2007 +0100
@@ -11,7 +11,7 @@
 #include <vector>
 
 #include "common.h"
-#include "kdtree.h"
+#include "container.h"
 #include "scene.h"
 
 using namespace std;
diff -r 28f6e8b9d5d1 -r fb170fccb19f src/Makefile
--- a/src/Makefile	Sun Dec 09 23:10:58 2007 +0100
+++ b/src/Makefile	Wed Dec 12 19:59:19 2007 +0100
@@ -6,7 +6,7 @@
 
 vpath %.cc $(ROOT)/src
 vpath %.h $(ROOT)/include
-LIBOBJS=raytracer.o scene.o noise.o container.o kdtree.o
+LIBOBJS=raytracer.o scene.o noise.o container.o kdtree.o octree.o
 CCFLAGS+=-I$(ROOT)/include
 
 ### Targets ###
@@ -45,6 +45,7 @@
 scene.o: scene.cc scene.h vector.h noise.h common.h
 container.o: container.cc container.h scene.h common.h
 kdtree.o: kdtree.cc kdtree.h scene.h common.h
+octree.o: octree.cc octree.h scene.h common.h
 raytracer.o: raytracer.cc raytracer.h scene.h vector.h noise.h common.h
 
 raytracermodule.o: raytracermodule.cc $(LIBOBJS)
diff -r 28f6e8b9d5d1 -r fb170fccb19f src/kdtree.cc
--- a/src/kdtree.cc	Sun Dec 09 23:10:58 2007 +0100
+++ b/src/kdtree.cc	Wed Dec 12 19:59:19 2007 +0100
@@ -379,7 +379,10 @@
 		for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)
 			if (*shape != origin_shape && (*shape)->intersect(ray, dist)
 			&& dist >= enPt->t)
+			{
 				nearest_shape = *shape;
+				nearest_distance = dist;
+			}
 
 		delete enPt;
 
@@ -390,7 +393,6 @@
 				delete st.back();
 				st.pop_back();
 			}
-			nearest_distance = dist;
 			return nearest_shape;
 		}
 
diff -r 28f6e8b9d5d1 -r fb170fccb19f src/octree.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/octree.cc	Wed Dec 12 19:59:19 2007 +0100
@@ -0,0 +1,325 @@
+/*
+ * Pyrit Ray Tracer
+ * file: octree.cc
+ *
+ * Radek Brich, 2006-2007
+ */
+
+#include "octree.h"
+
+OctreeNode::~OctreeNode()
+{
+	if (shapes != NULL)
+		delete shapes;
+	else
+		delete[] children;
+}
+
+void OctreeNode::subdivide(BBox bbox, int maxdepth)
+{
+	if (maxdepth <= 0 || shapes->size() <= 4)
+		return;
+
+	// make children
+	children = new OctreeNode[8];
+
+	// evaluate centres for axes
+	const Float xsplit = (bbox.L.x + bbox.H.x)*0.5;
+	const Float ysplit = (bbox.L.y + bbox.H.y)*0.5;
+	const Float zsplit = (bbox.L.z + bbox.H.z)*0.5;
+
+	// set bounding boxes for children
+	BBox childbb[8] = {bbox, bbox, bbox, bbox, bbox, bbox, bbox, bbox};
+	for (int i = 0; i < 4; i++)
+	{
+		// this is little obfuscated, so on right are listed affected children
+		// the idea is to cut every axis once per child, making 8 combinations
+		childbb[i].H.x = xsplit;	// 0,1,2,3
+		childbb[i+4].L.x = xsplit;  // 4,5,6,7
+		childbb[i+(i>>1<<1)].H.y = ysplit;	// 0,1,4,5
+		childbb[i+(i>>1<<1)+2].L.y = ysplit;// 2,3,6,7
+		childbb[i<<1].H.z = zsplit;		// 0,2,4,6
+		childbb[(i<<1)+1].L.z = zsplit; // 1,3,5,7
+	}
+
+	// distribute shapes to children
+	ShapeList::iterator sh;
+	BBox shbb;
+	int child, both;
+	unsigned int shapenum = 0;
+	for (sh = shapes->begin(); sh != shapes->end(); sh++)
+	{
+		child = 0;
+		both = 0;
+		shbb = (*sh)->get_bbox();
+
+		if (shbb.L.x >= xsplit)
+			child |= 4; //right
+		else if (shbb.H.x > xsplit)
+			both |= 4; // both
+		// for left, do nothing
+
+		if (shbb.L.y >= ysplit)
+			child |= 2;
+		else if (shbb.H.y > ysplit)
+			both |= 2;
+
+		if (shbb.L.z >= zsplit)
+			child |= 1;
+		else if (shbb.H.z > zsplit)
+			both |= 1;
+
+		if (!both)
+		{
+			getChild(child)->addShape(*sh);
+			shapenum++;
+		}
+		else
+		{
+			// shape goes to more than one child
+			if (both == 7)
+			{
+				for (int i = 0; i < 8; i++)
+					getChild(i)->addShape(*sh);
+				shapenum += 8;
+			}
+			else if (both == 3 || both >= 5)
+			{
+				if (both == 3)
+				{
+					for (int i = 0; i < 4; i++)
+						getChild(child + i)->addShape(*sh);
+				}
+				else if (both == 5)
+				{
+					for (int i = 0; i < 4; i++)
+						getChild(child + i+(i>>1<<1))->addShape(*sh);
+				}
+				else if (both == 6)
+				{
+					for (int i = 0; i < 4; i++)
+						getChild(child + (i<<1))->addShape(*sh);
+				}
+				shapenum += 4;
+			}
+			else
+			{
+				getChild(child)->addShape(*sh);
+				getChild(child+both)->addShape(*sh);
+				shapenum += 2;
+			}
+		}
+	}
+
+	if (shapes->size() <= 8 && shapenum > 2*shapes->size())
+	{
+		// bad subdivision, revert
+		delete[] children;
+		return;
+	}
+
+	// remove shapes and set this node to non-leaf
+	delete shapes;
+	shapes = NULL;
+
+	// recursive subdivision
+	for (int i = 0; i < 8; i++)
+		children[i].subdivide(childbb[i], maxdepth-1);
+}
+
+void Octree::build()
+{
+	dbgmsg(1, "* building octree\n");
+	root = new OctreeNode();
+	ShapeList::iterator shape;
+	for (shape = shapes.begin(); shape != shapes.end(); shape++)
+		root->addShape(*shape);
+
+	root->subdivide(bbox, max_depth);
+	built = true;
+}
+
+static inline int first_node(const Float tx0, const Float ty0, const Float tz0,
+	const Float txm, const Float tym, const Float tzm)
+{
+	int res = 0;
+	if (tx0 > ty0)
+	{
+		if (tx0 > tz0)
+		{ // YZ
+			if (tym < tx0)
+				res |= 2;
+			if (tzm < tx0)
+				res |= 1;
+		}
+		else
+		{ // XY
+			if (txm < tz0)
+				res |= 4;
+			if (tym < tz0)
+				res |= 2;
+		}
+	}
+	else
+	{
+		if (ty0 > tz0)
+		{ // XZ
+			if (txm < ty0)
+				res |= 4;
+			if (tzm < ty0)
+				res |= 1;
+			return res;
+		}
+		else
+		{ // XY
+			if (txm < tz0)
+				res |= 4;
+			if (tym < tz0)
+				res |= 2;
+		}
+	}
+	return res;
+}
+
+static inline int next_node(const Float txm, const int xnode,
+	const Float tym, const int ynode, const Float tzm, const int znode)
+{
+	if (txm < tym)
+	{
+		if (txm < tzm)
+			return xnode;
+		else
+			return znode;
+	}
+	else
+	{
+		if (tym < tzm)
+			return ynode;
+		else
+			return znode;
+	}
+}
+
+static Shape *proc_subtree(const int a, const Float tx0, const Float ty0, const Float tz0,
+	const Float tx1, const Float ty1, const Float tz1, OctreeNode *node,
+	const Shape *origin_shape, const Ray &ray, Float &nearest_distance)
+{
+	Float txm, tym, tzm;
+	int curr_node;
+
+	// if ray does not intersect this node
+	if (tx1 < 0.0 || ty1 < 0.0 || tz1 < 0.0)
+		return NULL;
+
+	if (node->isLeaf())
+	{
+		Shape *nearest_shape = NULL;
+		ShapeList::iterator shape;
+		Float mindist = max(max(tx0,ty0),tz0);
+		Float dist = min(min(min(tx1,ty1),tz1),nearest_distance);
+		for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)
+			if (*shape != origin_shape && (*shape)->intersect(ray, dist) && dist >= mindist)
+			{
+				nearest_shape = *shape;
+				nearest_distance = dist;
+			}
+		return nearest_shape;
+	}
+
+	txm = 0.5 * (tx0+tx1);
+	tym = 0.5 * (ty0+ty1);
+	tzm = 0.5 * (tz0+tz1);
+
+	curr_node = first_node(tx0,ty0,tz0,txm,tym,tzm);
+	Shape *shape = NULL;
+	while (curr_node < 8)
+	{
+		switch (curr_node)
+		{
+			case 0:
+				shape =proc_subtree (a,tx0,ty0,tz0,txm,tym,tzm,node->getChild(a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(txm, 4, tym, 2, tzm, 1);
+				break;
+			case 1:
+				shape =proc_subtree (a,tx0,ty0,tzm,txm,tym,tz1,node->getChild(1^a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(txm, 5, tym, 3, tz1, 8);
+				break;
+			case 2:
+				shape =proc_subtree (a,tx0,tym,tz0,txm,ty1,tzm,node->getChild(2^a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(txm, 6, ty1, 8, tzm, 3);
+				break;
+			case 3:
+				shape =proc_subtree (a,tx0,tym,tzm,txm,ty1,tz1,node->getChild(3^a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(txm, 7, ty1, 8, tz1, 8);
+				break;
+			case 4:
+				shape =proc_subtree (a,txm,ty0,tz0,tx1,tym,tzm,node->getChild(4^a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(tx1, 8, tym, 6, tzm, 5);
+				break;
+			case 5:
+				shape =proc_subtree (a,txm,ty0,tzm,tx1,tym,tz1,node->getChild(5^a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(tx1, 8, tym, 7, tz1, 8);
+				break;
+			case 6:
+				shape =proc_subtree (a,txm,tym,tz0,tx1,ty1,tzm,node->getChild(6^a), origin_shape, ray, nearest_distance);
+				curr_node = next_node(tx1, 8, ty1, 8, tzm, 7);
+				break;
+			case 7:
+				shape =proc_subtree (a,txm,tym,tzm,tx1,ty1,tz1,node->getChild(7^a), origin_shape, ray, nearest_distance);
+				curr_node = 8;
+				break;
+		}
+		if (shape != NULL)
+			return shape;
+	}
+	return NULL;
+}
+
+/*
+traversal algorithm paper as described in paper
+"An Efficient Parametric Algorithm for Octree Traversal"
+by J. Revelles, C. Urena and M. Lastra.
+*/
+Shape * Octree::nearest_intersection(const Shape *origin_shape, const Ray &ray,
+		Float &nearest_distance)
+{
+	/* if we have no tree, fall back to naive test */
+	if (!built)
+		return Container::nearest_intersection(origin_shape, ray, nearest_distance);
+
+	int a = 0;
+	Vector3 ro = ray.o;
+	Vector3 rdir = ray.dir;
+
+	if (rdir.x < 0.0)
+	{
+		ro.x = (bbox.L.x+bbox.H.x) - ro.x;
+		rdir.x = -rdir.x;
+		a |= 4;
+	}
+	if (rdir.y < 0.0)
+	{
+		ro.y = (bbox.L.y+bbox.H.y) - ro.y;
+		rdir.y = -rdir.y;
+		a |= 2;
+	}
+	if (rdir.z < 0.0)
+	{
+		ro.z = (bbox.L.z+bbox.H.z) - ro.z;
+		rdir.z = -rdir.z;
+		a |= 1;
+	}
+	Float tx0 = (bbox.L.x - ro.x) / rdir.x;
+	Float tx1 = (bbox.H.x - ro.x) / rdir.x;
+	Float ty0 = (bbox.L.y - ro.y) / rdir.y;
+	Float ty1 = (bbox.H.y - ro.y) / rdir.y;
+	Float tz0 = (bbox.L.z - ro.z) / rdir.z;
+	Float tz1 = (bbox.H.z - ro.z) / rdir.z;
+
+	//Octree *node = root;
+	if (max(max(tx0,ty0),tz0) < min (min(tx1,ty1),tz1))
+		return proc_subtree(a,tx0,ty0,tz0,tx1,ty1,tz1,root,
+			origin_shape, ray, nearest_distance);
+	else
+		return NULL;
+}
diff -r 28f6e8b9d5d1 -r fb170fccb19f src/raytracermodule.cc
--- a/src/raytracermodule.cc	Sun Dec 09 23:10:58 2007 +0100
+++ b/src/raytracermodule.cc	Wed Dec 12 19:59:19 2007 +0100
@@ -9,8 +9,9 @@
 #include <Python.h>
 
 #include <vector>
+#include "raytracer.h"
 #include "scene.h"
-#include "raytracer.h"
+#include "kdtree.h"
 
 //=========================== Light Source Object ===========================