# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1197500891 -3600
# Node ID b490093b0ac3a4de74f5b10faed7e63748eb41d4
# Parent  fb170fccb19f10f10712611e92a02ee436daa2a3
new virtual Shape::intersect_bbox
implementation of triangle-AABB intersection
octree building updated and simplified with help of this new method
octree made default for Python, it's currently much faster than kd-tree (both building and traversal)

diff -r fb170fccb19f -r b490093b0ac3 ccdemos/realtime_bunny.cc
--- a/ccdemos/realtime_bunny.cc	Wed Dec 12 19:59:19 2007 +0100
+++ b/ccdemos/realtime_bunny.cc	Thu Dec 13 00:08:11 2007 +0100
@@ -30,7 +30,6 @@
 	while (token != "end_header")
 	{
 		f >> token;
-		cout << token << endl;
 		if (token == "element")
 		{
 			f >> token;
@@ -146,7 +145,7 @@
 	render_buffer = (Float *) malloc(w*h*3*sizeof(Float));
 
 	rt.setThreads(2);
-	rt.setMaxDepth(3);
+	rt.setMaxDepth(0);
 
 	Octree top;
 	rt.setTop(&top);
@@ -160,13 +159,13 @@
 	//rt.addlight(&light2);
 
 	Material mat(Colour(0.9, 0.9, 0.9));
-	load_ply("../models/bunny/bun_zipper_res4.ply", &mat, 29);
+	load_ply("../models/bunny/bun_zipper.ply", &mat, 29);
 
 	rt.setCamera(&cam);
 	cam.setEye(Vector3(0,0,10));
 
 	/* build kd-tree */
-	top.setMaxDepth(3);
+	top.setMaxDepth(8);
 	top.optimize();
 
 	/* loop... */
diff -r fb170fccb19f -r b490093b0ac3 ccdemos/realtime_dragon.cc
--- a/ccdemos/realtime_dragon.cc	Wed Dec 12 19:59:19 2007 +0100
+++ b/ccdemos/realtime_dragon.cc	Thu Dec 13 00:08:11 2007 +0100
@@ -143,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));
diff -r fb170fccb19f -r b490093b0ac3 include/octree.h
--- a/include/octree.h	Wed Dec 12 19:59:19 2007 +0100
+++ b/include/octree.h	Thu Dec 13 00:08:11 2007 +0100
@@ -38,7 +38,7 @@
 	bool built;
 	int max_depth;
 public:
-	Octree() : Container(), root(NULL), built(false), max_depth(4) {};
+	Octree() : Container(), root(NULL), built(false), max_depth(10) {};
 	~Octree() { if (root) delete root; };
 	void addShape(Shape* aShape) { Container::addShape(aShape); built = false; };
 	Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
diff -r fb170fccb19f -r b490093b0ac3 include/scene.h
--- a/include/scene.h	Wed Dec 12 19:59:19 2007 +0100
+++ b/include/scene.h	Thu Dec 13 00:08:11 2007 +0100
@@ -130,6 +130,9 @@
 	// all intersections (only for CSG)
 	virtual bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const = 0;
 
+	// intersection with AABB
+	virtual bool intersect_bbox(const BBox &bbox) const = 0;
+
 	// normal at point P
 	virtual const Vector3 normal(const Vector3 &P) const = 0;
 
@@ -153,6 +156,7 @@
 		center(acenter), radius(aradius) { material = amaterial; }
 	bool intersect(const Ray &ray, Float &dist) const;
 	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const;
+	bool intersect_bbox(const BBox &bbox) const { return true; };
 	const Vector3 normal(const Vector3 &P) const { return (P - center) * inv_radius; };
 	BBox get_bbox() const;
 };
@@ -171,6 +175,7 @@
 	};
 	bool intersect(const Ray &ray, Float &dist) const;
 	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const { return false; };
+	bool intersect_bbox(const BBox &bbox) const { return true; };
 	const Vector3 normal(const Vector3 &P) const;
 	BBox get_bbox() const { return BBox(L, H); };
 };
@@ -235,6 +240,7 @@
 	Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial);
 	bool intersect(const Ray &ray, Float &dist) const;
 	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const {return false;};
+	bool intersect_bbox(const BBox &bbox) const;
 	const Vector3 normal(const Vector3 &P) const { return (smooth ? smooth_normal(P) : N); };
 	const Vector3 getNormal() const { return N; };
 	void setSmooth() { smooth = true; };//(typeid(*A) == typeid(*B) == typeid(*C) == typeid(NormalVertex)); };
diff -r fb170fccb19f -r b490093b0ac3 src/octree.cc
--- a/src/octree.cc	Wed Dec 12 19:59:19 2007 +0100
+++ b/src/octree.cc	Thu Dec 13 00:08:11 2007 +0100
@@ -44,71 +44,15 @@
 
 	// 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++)
+			if ((*sh)->intersect_bbox(childbb[i]))
 			{
-				for (int i = 0; i < 8; i++)
-					getChild(i)->addShape(*sh);
-				shapenum += 8;
+				getChild(i)->addShape(*sh);
+				shapenum++;
 			}
-			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())
@@ -316,7 +260,6 @@
 	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);
diff -r fb170fccb19f -r b490093b0ac3 src/raytracermodule.cc
--- a/src/raytracermodule.cc	Wed Dec 12 19:59:19 2007 +0100
+++ b/src/raytracermodule.cc	Thu Dec 13 00:08:11 2007 +0100
@@ -11,7 +11,7 @@
 #include <vector>
 #include "raytracer.h"
 #include "scene.h"
-#include "kdtree.h"
+#include "octree.h"
 
 //=========================== Light Source Object ===========================
 
@@ -571,7 +571,7 @@
 	v->raytracer = new Raytracer();
 	v->children = new vector<PyObject*>();
 	v->raytracer->setCamera(new Camera());
-	v->raytracer->setTop(new KdTree());
+	v->raytracer->setTop(new Octree());
 
 	return (PyObject*)v;
 }
diff -r fb170fccb19f -r b490093b0ac3 src/scene.cc
--- a/src/scene.cc	Wed Dec 12 19:59:19 2007 +0100
+++ b/src/scene.cc	Thu Dec 13 00:08:11 2007 +0100
@@ -329,6 +329,136 @@
 #endif
 }
 
+bool Triangle::intersect_bbox(const BBox &bbox) const
+{
+	const Vector3 boxcenter = (bbox.L+bbox.H)*0.5;
+	const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5;
+	const Vector3 v0 = A->P - boxcenter;
+	const Vector3 v1 = B->P - boxcenter;
+	const Vector3 v2 = C->P - boxcenter;
+	const Vector3 e0 = v1-v0;
+	const Vector3 e1 = v2-v1;
+	const Vector3 e2 = v0-v2;
+
+	Float fex = fabsf(e0.x);
+	Float fey = fabsf(e0.y);
+	Float fez = fabsf(e0.z);
+
+	Float p0,p1,p2,min,max,rad;
+
+	p0 = e0.z*v0.y - e0.y*v0.z;
+	p2 = e0.z*v2.y - e0.y*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = -e0.z*v0.x + e0.x*v0.z;
+	p2 = -e0.z*v2.x + e0.x*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p1 = e0.y*v1.x - e0.x*v1.y;
+	p2 = e0.y*v2.x - e0.x*v2.y;
+	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
+	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
+	if(min>rad || max<-rad) return false;
+
+	fex = fabsf(e1.x);
+	fey = fabsf(e1.y);
+	fez = fabsf(e1.z);
+
+	p0 = e1.z*v0.y - e1.y*v0.z;
+	p2 = e1.z*v2.y - e1.y*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = -e1.z*v0.x + e1.x*v0.z;
+	p2 = -e1.z*v2.x + e1.x*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = e1.y*v0.x - e1.x*v0.y;
+	p1 = e1.y*v1.x - e1.x*v1.y;
+	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
+	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
+	if(min>rad || max<-rad) return false;
+
+	fex = fabsf(e2.x);
+	fey = fabsf(e2.y);
+	fez = fabsf(e2.z);
+
+	p0 = e2.z*v0.y - e2.y*v0.z;
+	p1 = e2.z*v1.y - e2.y*v1.z;
+	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
+	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = -e2.z*v0.x + e2.x*v0.z;
+	p1 = -e2.z*v1.x + e2.x*v1.z;
+	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
+	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p1 = e2.y*v1.x - e2.x*v1.y;
+	p2 = e2.y*v2.x - e2.x*v2.y;
+	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
+	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
+	if(min>rad || max<-rad) return false;
+
+	/* test overlap in the {x,y,z}-directions */
+	/* test in X-direction */
+	min = v0.x;
+	if (v1.x < min) min = v1.x;
+	if (v2.x < min) min = v2.x;
+	max = v0.x;
+	if (v1.x > max) max = v1.x;
+	if (v2.x > max) max = v2.x;
+	if(min>boxhalfsize.x || max<-boxhalfsize.x) return false;
+
+	/* test in Y-direction */
+	min = v0.y;
+	if (v1.y < min) min = v1.y;
+	if (v2.y < min) min = v2.y;
+	max = v0.y;
+	if (v1.y > max) max = v1.y;
+	if (v2.y > max) max = v2.y;
+	if(min>boxhalfsize.y || max<-boxhalfsize.y) return false;
+
+	/* test in Z-direction */
+	min = v0.z;
+	if (v1.z < min) min = v1.z;
+	if (v2.z < min) min = v2.z;
+	max = v0.z;
+	if (v1.z > max) max = v1.z;
+	if (v2.z > max) max = v2.z;
+	if(min>boxhalfsize.z || max<-boxhalfsize.z) return false;
+
+	/*  test if the box intersects the plane of the triangle */
+	Vector3 vmin,vmax;
+	Float v;
+	for(int q=0;q<3;q++)
+	{
+		v=v0[q];
+		if(N[q]>0.0f)
+		{
+			vmin.cell[q]=-boxhalfsize[q] - v;
+			vmax.cell[q]= boxhalfsize[q] - v;
+		}
+		else
+		{
+			vmin.cell[q]= boxhalfsize[q] - v;
+			vmax.cell[q]=-boxhalfsize[q] - v;
+		}
+	}
+	if(dot(N,vmin)>0.0f) return false;
+	if(dot(N,vmax)>=0.0f) return true;
+
+	return false;
+}
+
 BBox Triangle::get_bbox() const
 {
 	BBox bbox = BBox();