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)
/* * 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; unsigned int shapenum = 0; for (sh = shapes->begin(); sh != shapes->end(); sh++) { for (int i = 0; i < 8; i++) if ((*sh)->intersect_bbox(childbb[i])) { getChild(i)->addShape(*sh); shapenum++; } } 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; 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;}