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)
/* * Pyrit Ray Tracer * file: kdtree.cc * * Radek Brich, 2006-2007 */#include <algorithm>#include <stack>#include "common.h"#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 traversalclass 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) {};};// ----------------------------------------KdNode::~KdNode(){ if (isLeaf()) delete shapes; else delete[] children;}void KdNode::subdivide(BBox bbox, int maxdepth){ if (maxdepth <= 0 || shapes->size() <= 2) { setLeaf(); return; } // choose split axis axis = 0; if (bbox.h() > bbox.w() && bbox.h() > bbox.d()) axis = 1; if (bbox.d() > bbox.w() && bbox.d() > bbox.h()) axis = 2; // *** find optimal split position SortableShapeList sslist(shapes, axis); sort(sslist.begin(), sslist.end()); 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.H.cell[axis])); } sort(splitlist.begin(), splitlist.end()); // remove duplicate splits splitlist.resize(unique(splitlist.begin(), splitlist.end()) - splitlist.begin()); // find all posible splits SplitList::iterator spl; int rest; for (spl = splitlist.begin(); spl != splitlist.end(); spl++) { for (sh = sslist.begin(), rest = sslist.size(); sh != sslist.end(); sh++, rest--) { if (sh->hasMark()) { spl->lnum++; continue; } // if shape is completely contained in split plane if (spl->pos == sh->bbox.L.cell[axis] == sh->bbox.H.cell[axis]) { if (spl->pos - bbox.L.cell[axis] < bbox.H.cell[axis] - spl->pos) { // left subcell is smaller -> if not empty, put shape here if (spl->lnum) spl->lnum++; else spl->rnum++; } else { // right subcell is smaller if (spl->rnum) spl->rnum++; else spl->lnum++; } sh->setMark(); } else // if shape is on left side of split plane if (sh->bbox.H.cell[axis] <= spl->pos) { spl->lnum++; sh->setMark(); } else // if shape occupies both sides of split plane if (sh->bbox.L.cell[axis] < spl->pos && sh->bbox.H.cell[axis] > spl->pos) { spl->lnum++; spl->rnum++; } else // if shape is on right side of split plane if (sh->bbox.L.cell[axis] >= spl->pos) { spl->rnum += rest; break; } } } // choose best split pos const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection Float SAV = 2*(bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost bool leaf = true; BBox lbb = bbox; BBox rbb = bbox; for (spl = splitlist.begin(); spl != splitlist.end(); spl++) { // calculate SAH cost of this split lbb.H.cell[axis] = spl->pos; rbb.L.cell[axis] = spl->pos; Float SAL = 2*(lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d()); Float SAR = 2*(rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d()); Float splitcost = K + SAL/SAV*(K+spl->lnum) + SAR/SAV*(K+spl->rnum); if (splitcost < cost) { leaf = false; cost = splitcost; split = spl->pos; } } if (leaf) { setLeaf(); return; }#if 0// export kd-tree as .obj for visualization// this would be hard to reconstruct later static ofstream F("kdtree.obj"); Vector3 v; static int f=0; v.cell[axis] = split; v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3]; v.cell[(axis+2)%3] = bbox.L.cell[(axis+2)%3]; F << "v " << v.x << " " << v.y << " " << v.z << endl; v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3]; v.cell[(axis+2)%3] = bbox.H.cell[(axis+2)%3]; F << "v " << v.x << " " << v.y << " " << v.z << endl; v.cell[(axis+1)%3] = bbox.H.cell[(axis+1)%3]; v.cell[(axis+2)%3] = bbox.H.cell[(axis+2)%3]; F << "v " << v.x << " " << v.y << " " << v.z << endl; v.cell[(axis+1)%3] = bbox.H.cell[(axis+1)%3]; v.cell[(axis+2)%3] = bbox.L.cell[(axis+2)%3]; F << "v " << v.x << " " << v.y << " " << v.z << endl; F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl; f += 4;#endif // split this node delete shapes; children = new KdNode[2]; int state = 0; for (sh = sslist.begin(); sh != sslist.end(); sh++) { // if shape is on left side of split plane if (state == 1) { // only right children[1].addShape(sh->shape); continue; } if (state == 0) { if (sh->bbox.H.cell[axis] < split) { // left children[0].addShape(sh->shape); } else if (sh->bbox.H.cell[axis] > split) { if (sh->bbox.L.cell[axis] < split) { // both children[0].addShape(sh->shape); children[1].addShape(sh->shape); } else { // right children[1].addShape(sh->shape); state = 1; } } else { // H == split if (sh->bbox.L.cell[axis] < split) { // left children[0].addShape(sh->shape); } else { // contained if (split - bbox.L.cell[axis] < bbox.H.cell[axis] - split) { // left subcell is smaller -> if not empty, put shape here if (children[0].shapes->size()) children[0].addShape(sh->shape); else children[1].addShape(sh->shape); } else { // right subcell is smaller if (children[1].shapes->size()) children[1].addShape(sh->shape); else children[0].addShape(sh->shape); } } } } } lbb.H.cell[axis] = split; rbb.L.cell[axis] = split; children[0].subdivide(lbb, maxdepth-1); children[1].subdivide(rbb, maxdepth-1);}void KdTree::build(){ dbgmsg(1, "* building kd-tree\n"); root = new KdNode(); ShapeList::iterator shape; for (shape = shapes.begin(); shape != shapes.end(); shape++) root->addShape(*shape); root->subdivide(bbox, max_depth); 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; /* pointers to the far child node and current node */ KdNode *farchild, *node; node = root; /* start from the kd-tree root node */ /* std vector is much faster than stack */ vector<StackElem*> st; 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 */ /* setup initial exit point in the stack */ StackElem *exPt = new StackElem(NULL, b, ray.o + ray.dir * b); st.push_back(exPt); /* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */ 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(); 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.o.cell[axis]) / ray.dir.cell[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 */ /* 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; nearest_distance = dist; } delete enPt; if (nearest_shape) { while (!st.empty()) { delete st.back(); st.pop_back(); } return nearest_shape; } enPt = exPt; /* retrieve the pointer to the next node, it is possible that ray traversal terminates */ node = enPt->node; st.pop_back(); } /* while */ delete enPt; /* ray leaves the scene */ return NULL;}// this should save whole kd-tree with triangles distributed into leavesvoid KdTree::save(ostream &str, KdNode *node){ if (!built) return; if (node == NULL) node = root; if (node->isLeaf()) str << "(leaf: " << node->shapes->size() << " shapes)"; else { str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L="; save(str, node->getLeftChild()); str << "; R="; save(str, node->getRightChild()); str << ";)"; }}// load kd-tree from file/streamvoid KdTree::load(istream &str, KdNode *node){}