--- a/src/kdtree.cc Sun Apr 20 19:27:59 2008 +0200
+++ b/src/kdtree.cc Mon Apr 21 08:47:36 2008 +0200
@@ -48,14 +48,12 @@
};
// stack element for kd-tree traversal
-class StackElem
+struct StackElem
{
-public:
KdNode* node; /* pointer to far child */
Float t; /* the entry/exit signed distance */
Vector3 pb; /* the coordinates of entry/exit point */
- StackElem(KdNode *anode, const Float &at, const Vector3 &apb):
- node(anode), t(at), pb(apb) {};
+ int prev;
};
// ----------------------------------------
@@ -63,15 +61,15 @@
KdNode::~KdNode()
{
if (isLeaf())
- delete shapes;
+ delete getShapes();
else
- delete[] children;
+ delete[] getLeftChild();
}
// kd-tree recursive build algorithm, inspired by PBRT (www.pbrt.org)
void KdNode::subdivide(BBox bounds, int maxdepth)
{
- if (maxdepth <= 0 || shapes->size() <= 2)
+ if (maxdepth <= 0 || getShapes()->size() <= 2)
{
setLeaf();
return;
@@ -87,7 +85,7 @@
// create sorted list of shape bounds (= find all posible splits)
vector<ShapeBound> edges[3];
ShapeList::iterator shape;
- for (shape = shapes->begin(); shape != shapes->end(); shape++)
+ for (shape = getShapes()->begin(); shape != getShapes()->end(); shape++)
{
BBox shapebounds = (*shape)->get_bbox();
for (int ax = 0; ax < 3; ax++)
@@ -103,12 +101,13 @@
const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
Float SAV = (bounds.w()*bounds.h() + // surface area of node
bounds.w()*bounds.d() + bounds.h()*bounds.d());
- Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost
+ Float cost = SAV * (K + getShapes()->size()); // initial cost = non-split cost
vector<ShapeBound>::iterator edge, splitedge = edges[2].end();
+ int axis = 0;
for (int ax = 0; ax < 3; ax++)
{
- int lnum = 0, rnum = shapes->size();
+ int lnum = 0, rnum = getShapes()->size();
BBox lbb = bounds;
BBox rbb = bounds;
for (edge = edges[ax].begin(); edge != edges[ax].end(); edge++)
@@ -166,22 +165,24 @@
#endif
// split this node
- delete shapes;
+ delete getShapes();
+
BBox lbb = bounds;
BBox rbb = bounds;
lbb.H.cell[axis] = split;
rbb.L.cell[axis] = split;
- children = new KdNode[2];
+ setChildren(new KdNode[2]);
+ setAxis(axis);
for (edge = edges[axis].begin(); edge != splitedge; edge++)
if (!edge->end && edge->shape->intersect_bbox(lbb))
- children[0].addShape(edge->shape);
+ getLeftChild()->addShape(edge->shape);
for (edge = splitedge; edge < edges[axis].end(); edge++)
if (edge->end && edge->shape->intersect_bbox(rbb))
- children[1].addShape(edge->shape);
+ getRightChild()->addShape(edge->shape);
- children[0].subdivide(lbb, maxdepth-1);
- children[1].subdivide(rbb, maxdepth-1);
+ getLeftChild()->subdivide(lbb, maxdepth-1);
+ getRightChild()->subdivide(rbb, maxdepth-1);
}
void KdTree::build()
@@ -211,40 +212,46 @@
/* pointers to the far child node and current node */
KdNode *farchild, *node;
- node = root; /* start from the kd-tree root node */
+ node = root;
- /* std vector is much faster than stack */
- vector<StackElem*> st;
+ StackElem stack[max_depth];
- StackElem *enPt = new StackElem(NULL, a,
- /* distinguish between internal and external origin of a ray*/
- a >= 0.0 ?
- ray.o + ray.dir * a : /* external */
- ray.o); /* internal */
+ int entry = 0, exit = 1;
+ stack[entry].t = a;
+
+ /* distinguish between internal and external origin of a ray*/
+ if (a >= 0.0)
+ stack[entry].pb = ray.o + ray.dir * a; /* external */
+ else
+ stack[entry].pb = ray.o; /* internal */
/* setup initial exit point in the stack */
- StackElem *exPt = new StackElem(NULL, b, ray.o + ray.dir * b);
- st.push_back(exPt);
+ stack[exit].t = b;
+ stack[exit].pb = ray.o + ray.dir * b;
+ stack[exit].node = NULL;
/* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */
+ Float splitVal;
+ int axis;
+ static const int mod3[] = {0,1,2,0,1};
+ const Vector3 invdir = 1 / ray.dir;
while (node)
{
- exPt = st.back();
/* loop until a leaf is found */
while (!node->isLeaf())
{
/* retrieve position of splitting plane */
- Float splitVal = node->getSplit();
- short axis = node->getAxis();
+ splitVal = node->getSplit();
+ axis = node->getAxis();
- if (enPt->pb[axis] <= splitVal)
+ if (stack[entry].pb[axis] <= splitVal)
{
- if (exPt->pb[axis] <= splitVal)
+ if (stack[exit].pb[axis] <= splitVal)
{ /* case N1, N2, N3, P5, Z2, and Z3 */
node = node->getLeftChild();
continue;
}
- if (exPt->pb[axis] == splitVal)
+ if (stack[exit].pb[axis] == splitVal)
{ /* case Z1 */
node = node->getRightChild();
continue;
@@ -254,14 +261,14 @@
node = node->getLeftChild();
}
else
- { /* (enPt->pb[axis] > splitVal) */
- if (splitVal < exPt->pb[axis])
+ { /* (stack[entry].pb[axis] > splitVal) */
+ if (splitVal < stack[exit].pb[axis])
{
/* case P1, P2, P3, and N5 */
node = node->getRightChild();
continue;
}
- /* case P4*/
+ /* case P4 */
farchild = node->getLeftChild();
node = node->getRightChild();
}
@@ -269,50 +276,48 @@
/* case P4 or N4 . . . traverse both children */
/* signed distance to the splitting plane */
- t = (splitVal - ray.o.cell[axis]) / ray.dir.cell[axis];
+ t = (splitVal - ray.o[axis]) * invdir[axis];
/* setup the new exit point and push it onto stack */
- exPt = new StackElem(farchild, t, Vector3());
- exPt->pb.cell[axis] = splitVal;
- exPt->pb.cell[(axis+1)%3] = ray.o.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
- exPt->pb.cell[(axis+2)%3] = ray.o.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
- st.push_back(exPt);
- } /* while */
+ register int tmp = exit;
+
+ exit++;
+ if (exit == entry)
+ exit++;
+ assert(exit < max_depth);
+
+ stack[exit].prev = tmp;
+ stack[exit].t = t;
+ stack[exit].node = farchild;
+ stack[exit].pb.cell[axis] = splitVal;
+ stack[exit].pb.cell[mod3[axis+1]] = ray.o.cell[mod3[axis+1]]
+ + t * ray.dir.cell[mod3[axis+1]];
+ stack[exit].pb.cell[mod3[axis+2]] = ray.o.cell[mod3[axis+2]]
+ + t * ray.dir.cell[mod3[axis+2]];
+ }
/* 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;
+ Float dist = stack[exit].t;
ShapeList::iterator shape;
- for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)
+ for (shape = node->getShapes()->begin(); shape != node->getShapes()->end(); shape++)
if (*shape != origin_shape && (*shape)->intersect(ray, dist)
- && dist >= enPt->t)
+ && dist >= stack[entry].t)
{
nearest_shape = *shape;
nearest_distance = dist;
}
- delete enPt;
-
if (nearest_shape)
- {
- while (!st.empty())
- {
- delete st.back();
- st.pop_back();
- }
return nearest_shape;
- }
- enPt = exPt;
+ entry = exit;
- /* retrieve the pointer to the next node, it is possible that ray traversal terminates */
- node = enPt->node;
- st.pop_back();
- } /* while */
-
- delete enPt;
+ /* retrieve the pointer to the next node,
+ it is possible that ray traversal terminates */
+ node = stack[entry].node;
+ exit = stack[entry].prev;
+ }
/* ray leaves the scene */
return NULL;
@@ -326,7 +331,7 @@
if (node == NULL)
node = root;
if (node->isLeaf())
- str << "(leaf: " << node->shapes->size() << " shapes)";
+ str << "(leaf: " << node->getShapes()->size() << " shapes)";
else
{
str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L=";