# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1195937741 -3600
# Node ID f4fcabf05785051a13f61362bf5ec3b00fa046de
# Parent  4d192e13ee848d23f0d412b8f1dea6fb3a27e9c4
kd-tree: traversal algorithm (KdTree::nearest_intersection)

diff -r 4d192e13ee84 -r f4fcabf05785 TODO
--- a/TODO	Fri Nov 23 16:14:38 2007 +0100
+++ b/TODO	Sat Nov 24 21:55:41 2007 +0100
@@ -4,6 +4,9 @@
  * transforms, camera
  * update Python binding, new classes
  * C++ demos
+ * rename:
+   - Ray.a -> Ray.o or Ray.origin
+   - BBox.R -> H
 
 New Classes?
 ============
diff -r 4d192e13ee84 -r f4fcabf05785 src/kdtree.cc
--- a/src/kdtree.cc	Fri Nov 23 16:14:38 2007 +0100
+++ b/src/kdtree.cc	Sat Nov 24 21:55:41 2007 +0100
@@ -1,7 +1,63 @@
 #include <algorithm>
+#include <stack>
 
 #include "kdtree.h"
 
+class SortableShape
+{
+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)
+	{
+		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(): 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
+struct StackElem
+{
+	KdNode* node; /* pointer to far child */
+	float t; /* the entry/exit signed distance */
+	Vector3 pb; /* the coordinates of entry/exit point */
+};
+
+// ----------------------------------------
+
 void Container::addShape(Shape* aShape)
 {
 	shapes.push_back(aShape);
@@ -238,6 +294,132 @@
 	built = true;
 }
 
+/* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */
+Shape *KdTree::nearest_intersection(const Shape *origin_shape, const Ray &ray,
+	float &nearest_distance)
+{
+	float a, b; /* entry/exit signed distance */
+	float t;    /* signed distance to the splitting plane */
+
+	/* if we have no tree, fall back to naive test */
+	if (!built)
+		return Container::nearest_intersection(origin_shape, ray, nearest_distance);
+
+	if (!bbox.intersect(ray, a, b))
+		return NULL;
+
+	stack<StackElem*> st;
+
+	/* pointers to the far child node and current node */
+	KdNode *farchild, *node;
+	node = root; /* start from the kd-tree root node */
+
+	StackElem *enPt = new StackElem();
+	enPt->t = a; /* set the signed distance */
+	enPt->node = NULL;
+
+	/* distinguish between internal and external origin */
+	if (a >= 0.0) /* a ray with external origin */
+		enPt->pb = ray.a + ray.dir * a;
+	else /* a ray with internal origin */
+		enPt->pb = ray.a;
+
+	/* setup initial exit point in the stack */
+	StackElem *exPt = new StackElem();
+	exPt->t = b;
+	exPt->pb = ray.a + ray.dir * b;
+	exPt->node = NULL; /* set termination flag */
+
+	st.push(enPt);
+
+	/* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */
+	while (node)
+	{
+		/* loop until a leaf is found */
+		while (!node->isLeaf())
+		{
+			/* retrieve position of splitting plane */
+			float splitVal = node->getSplit();
+			short axis = node->getAxis();
+
+			if (enPt->pb[axis] <= splitVal)
+			{
+				if (exPt->pb[axis] <= splitVal)
+				{ /* case N1, N2, N3, P5, Z2, and Z3 */
+					node = node->getLeftChild();
+					continue;
+				}
+				if (exPt->pb[axis] == splitVal)
+				{ /* case Z1 */
+					node = node->getRightChild();
+					continue;
+				}
+				/* case N4 */
+				farchild = node->getRightChild();
+				node = node->getLeftChild();
+			}
+			else
+			{ /* (enPt->pb[axis] > splitVal) */
+				if (splitVal < exPt->pb[axis])
+				{
+					/* case P1, P2, P3, and N5 */
+					node = node->getRightChild();
+					continue;
+				}
+				/* case P4*/
+				farchild = node->getLeftChild();
+				node = node->getRightChild();
+			}
+
+			/* case P4 or N4 . . . traverse both children */
+
+			/* signed distance to the splitting plane */
+			t = (splitVal - ray.a.cell[axis]) / ray.dir.cell[axis];
+
+			/* setup the new exit point */
+			st.push(exPt);
+			exPt = new StackElem();
+
+			/* push values onto the stack */
+			exPt->t = t;
+			exPt->node = farchild;
+			exPt->pb[axis] = splitVal;
+			exPt->pb[(axis+1)%3] = ray.a.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
+			exPt->pb[(axis+2)%3] = ray.a.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
+		} /* while */
+
+		/* current node is the leaf . . . empty or full */
+		/* "intersect ray with each object in the object list, discarding "
+		"those lying before stack[enPt].t or farther than stack[exPt].t" */
+		Shape *nearest_shape = NULL;
+		float dist = exPt->t;
+		ShapeList::iterator shape;
+		for (shape = node->shapes.begin(); shape != node->shapes.end(); shape++)
+			if (*shape != origin_shape && (*shape)->intersect(ray, dist)
+			&& dist >= enPt->t)
+				nearest_shape = *shape;
+
+		if (nearest_shape)
+		{
+			nearest_distance = dist;
+			return nearest_shape;
+		}
+
+		/* pop from the stack */
+		enPt = exPt; /* the signed distance intervals are adjacent */
+
+		/* retrieve the pointer to the next node, it is possible that ray traversal terminates */
+		node = enPt->node;
+
+		// pop
+		exPt = st.top();
+		st.pop();
+	} /* while */
+
+	/* ray leaves the scene */
+	return NULL;
+}
+
 // this should save whole kd-tree with triangles distributed into leaves
 void KdTree::save(ostream &str, KdNode *node)
 {
@@ -260,5 +442,5 @@
 // load kd-tree from file/stream
 void KdTree::load(istream &str, KdNode *node)
 {
-	
+
 }
diff -r 4d192e13ee84 -r f4fcabf05785 src/kdtree.h
--- a/src/kdtree.h	Fri Nov 23 16:14:38 2007 +0100
+++ b/src/kdtree.h	Sat Nov 24 21:55:41 2007 +0100
@@ -13,51 +13,6 @@
 {
 };
 
-class SortableShape
-{
-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)
-	{
-		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(): 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>
-{
-};
-
 class Container
 {
 protected:
@@ -68,7 +23,7 @@
 	virtual ~Container() {};
 	virtual void addShape(Shape* aShape);
 	//void addShapeNoExtend(Shape* aShape) { shapes.push_back(aShape); };
-        virtual Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
+	virtual Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
 		float &nearest_distance);
 	virtual void optimize() {};
 };
@@ -108,8 +63,10 @@
 public:
 	KdTree() : Container(), built(false) {};
 	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 optimize() { build(); };
 	void save(ostream &str, KdNode *node = NULL);
 	void load(istream &str, KdNode *node = NULL);
 };
diff -r 4d192e13ee84 -r f4fcabf05785 src/raytracer.cc
--- a/src/raytracer.cc	Fri Nov 23 16:14:38 2007 +0100
+++ b/src/raytracer.cc	Sat Nov 24 21:55:41 2007 +0100
@@ -212,7 +212,7 @@
 
 	printf("* building kd-tree\n");
 	//cout << endl;
-	static_cast<KdTree*>(top)->optimize();
+	top->optimize();
 	//static_cast<KdTree*>(top)->save(cout);
 	//cout << endl;
 
diff -r 4d192e13ee84 -r f4fcabf05785 src/scene.cc
--- a/src/scene.cc	Fri Nov 23 16:14:38 2007 +0100
+++ b/src/scene.cc	Sat Nov 24 21:55:41 2007 +0100
@@ -6,8 +6,48 @@
  */
 
 #include <math.h>
+#include <float.h>
+
 #include "scene.h"
 
+/* http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm */
+bool BBox::intersect(const Ray &ray, float &a, float &b)
+{
+	float tnear = -FLT_MAX;
+	float tfar = FLT_MAX;
+	float t1, t2;
+
+	for (int i = 0; i < 3; i++)
+	{
+		if (ray.dir.cell[i] == 0) {
+			/* ray is parallel to these planes */
+			if (ray.a.cell[i] < L.cell[i] || ray.a.cell[i] > R.cell[i])
+				return false;
+		} else
+		{
+			/* compute the intersection distance of the planes */
+			t1 = (L.cell[i] - ray.a.cell[i]) / ray.dir.cell[i];
+			t2 = (R.cell[i] - ray.a.cell[i]) / ray.dir.cell[i];
+
+			if (t1 > t2)
+				swap(t1, t2);
+
+			if (t1 > tnear)
+				tnear = t1; /* want largest Tnear */
+			if (t2 < tfar)
+				tfar = t2; /* want smallest Tfar */
+			if (tnear > tfar)
+				return false; /* box missed */
+			if (tfar < 0)
+				return false; /* box is behind ray */
+		}
+	}
+
+	a = tnear;
+	b = tfar;
+	return true;
+}
+
 bool Sphere::intersect(const Ray &ray, float &dist)
 {
 	Vector3 V = ((Ray)ray).a - center;
@@ -72,7 +112,7 @@
 }
 
 BBox Sphere::get_bbox()
-{	
+{
 	BBox bbox = BBox();
 	bbox.L.x = center.x - radius;
 	bbox.L.y = center.y - radius;
@@ -98,7 +138,7 @@
 }
 
 BBox Plane::get_bbox()
-{	
+{
 	return BBox();
 }
 
diff -r 4d192e13ee84 -r f4fcabf05785 src/scene.h
--- a/src/scene.h	Fri Nov 23 16:14:38 2007 +0100
+++ b/src/scene.h	Sat Nov 24 21:55:41 2007 +0100
@@ -16,6 +16,14 @@
 
 using namespace std;
 
+class Ray
+{
+public:
+	Vector3 a, dir;
+	Ray(const Vector3 &aa, const Vector3 &adir):
+		a(aa), dir(adir) {};
+};
+
 /* axis-aligned bounding box */
 class BBox
 {
@@ -26,14 +34,7 @@
 	float w() { return R.x-L.x; };
 	float h() { return R.y-L.y; };
 	float d() { return R.z-L.z; };
-};
-
-class Ray
-{
-public:
-	Vector3 a, dir;
-	Ray(const Vector3 &aa, const Vector3 &adir):
-		a(aa), dir(adir) {};
+	bool intersect(const Ray &ray, float &a, float &b);
 };
 
 class Light
diff -r 4d192e13ee84 -r f4fcabf05785 src/vector.h
--- a/src/vector.h	Fri Nov 23 16:14:38 2007 +0100
+++ b/src/vector.h	Sat Nov 24 21:55:41 2007 +0100
@@ -94,9 +94,9 @@
 	};
 
 	// product of vector and scalar
-	Vector3 operator*(const float &f)
+	friend Vector3 operator*(const Vector3 &v, const float &f)
 	{
-		return Vector3(f * x, f * y, f * z);
+		return Vector3(f * v.x, f * v.y, f * v.z);
 	}
 
 	friend Vector3 operator*(const float &f, Vector3 &v)