diff -r b490093b0ac3 -r 5f954c0d34fc src/octree.cc --- a/src/octree.cc Thu Dec 13 00:08:11 2007 +0100 +++ b/src/octree.cc Fri Dec 14 00:05:54 2007 +0100 @@ -83,50 +83,32 @@ 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) + +/******************************************************* +octree traversal algorithm as described in paper +"An Efficient Parametric Algorithm for Octree Traversal" +by J. Revelles, C. Urena and M. Lastra. + +see revision 37 for original recursive version +*******************************************************/ + +struct OctreeTravState { - 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; -} + Float tx0,ty0,tz0,tx1,ty1,tz1,txm,tym,tzm; + OctreeNode *node; + int next; + OctreeTravState() {}; + OctreeTravState( + const Float atx0, const Float aty0, const Float atz0, + const Float atx1, const Float aty1, const Float atz1, + const Float atxm, const Float atym, const Float atzm, + OctreeNode *const anode, const int anext): + tx0(atx0), ty0(aty0), tz0(atz0), tx1(atx1), ty1(aty1), tz1(atz1), + txm(atxm), tym(atym), tzm(atzm), node(anode), next(anext) {}; +}; -static inline int next_node(const Float txm, const int xnode, - const Float tym, const int ynode, const Float tzm, const int znode) +inline const int &next_node(const Float &txm, const int &xnode, + const Float &tym, const int &ynode, const Float &tzm, const int &znode) { if (txm < tym) { @@ -144,86 +126,6 @@ } } -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) { @@ -231,9 +133,23 @@ if (!built) return Container::nearest_intersection(origin_shape, ray, nearest_distance); + OctreeTravState st[max_depth+1]; + register OctreeTravState *st_cur = st; + +# define node st_cur->node +# define tx0 st_cur->tx0 +# define ty0 st_cur->ty0 +# define tz0 st_cur->tz0 +# define tx1 st_cur->tx1 +# define ty1 st_cur->ty1 +# define tz1 st_cur->tz1 +# define txm st_cur->txm +# define tym st_cur->tym +# define tzm st_cur->tzm + int a = 0; - Vector3 ro = ray.o; - Vector3 rdir = ray.dir; + Vector3 ro(ray.o); + Vector3 rdir(1.0/ray.dir.x, 1.0/ray.dir.y, 1.0/ray.dir.z); if (rdir.x < 0.0) { @@ -253,16 +169,163 @@ 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; + + tx0 = (bbox.L.x - ro.x) * rdir.x; + tx1 = (bbox.H.x - ro.x) * rdir.x; + ty0 = (bbox.L.y - ro.y) * rdir.y; + ty1 = (bbox.H.y - ro.y) * rdir.y; + tz0 = (bbox.L.z - ro.z) * rdir.z; + tz1 = (bbox.H.z - ro.z) * rdir.z; + + if (max3(tx0,ty0,tz0) > min3(tx1,ty1,tz1)) + return NULL; + + node = root; + st_cur->next = -1; + + Shape *nearest_shape = NULL; + while (nearest_shape == NULL) + { + if (st_cur->next == -1) + { + st_cur->next = 8; + + // if ray does intersect this node + if (!(tx1 < 0.0 || ty1 < 0.0 || tz1 < 0.0)) + { + if (node->isLeaf()) + { + ShapeList::iterator shape; + register Float mindist = max3(tx0,ty0,tz0); + /* correct & slow */ + //Float dist = min(min3(tx1,ty1,tz1),nearest_distance); + /* faster */ + register Float dist = 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; + } + } + else + { + txm = 0.5 * (tx0+tx1); + tym = 0.5 * (ty0+ty1); + tzm = 0.5 * (tz0+tz1); - 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; + // first node + st_cur->next = 0; + if (tx0 > ty0) + { + if (tx0 > tz0) + { // YZ + if (tym < tx0) + st_cur->next |= 2; + if (tzm < tx0) + st_cur->next |= 1; + } + else + { // XY + if (txm < tz0) + st_cur->next |= 4; + if (tym < tz0) + st_cur->next |= 2; + } + } + else + { + if (ty0 > tz0) + { // XZ + if (txm < ty0) + st_cur->next |= 4; + if (tzm < ty0) + st_cur->next |= 1; + } + else + { // XY + if (txm < tz0) + st_cur->next |= 4; + if (tym < tz0) + st_cur->next |= 2; + } + } + } + } + } + + while (st_cur->next == 8) + { + // pop state from stack + if (st_cur == st) + return NULL; // nothing to pop, finish + --st_cur; + } + + // push current state + *(st_cur+1) = *st_cur; + ++st_cur; + + switch (st_cur->next) + { + case 0: + tx1 = txm; + ty1 = tym; + tz1 = tzm; + node = node->getChild(a); + (st_cur-1)->next = next_node(txm, 4, tym, 2, tzm, 1); + break; + case 1: + tz0 = tzm; + tx1 = txm; + ty1 = tym; + node = node->getChild(1^a); + (st_cur-1)->next = next_node(txm, 5, tym, 3, tz1, 8); + break; + case 2: + ty0 = tym; + tx1 = txm; + tz1 = tzm; + node = node->getChild(2^a); + (st_cur-1)->next = next_node(txm, 6, ty1, 8, tzm, 3); + break; + case 3: + ty0 = tym; + tz0 = tzm; + tx1 = txm; + node = node->getChild(3^a); + (st_cur-1)->next = next_node(txm, 7, ty1, 8, tz1, 8); + break; + case 4: + tx0 = txm; + ty1 = tym; + tz1 = tzm; + node = node->getChild(4^a); + (st_cur-1)->next = next_node(tx1, 8, tym, 6, tzm, 5); + break; + case 5: + tx0 = txm; + tz0 = tzm; + ty1 = tym; + node = node->getChild(5^a); + (st_cur-1)->next = next_node(tx1, 8, tym, 7, tz1, 8); + break; + case 6: + tx0 = txm; + ty0 = tym; + tz1 = tzm; + node = node->getChild(6^a); + (st_cur-1)->next = next_node(tx1, 8, ty1, 8, tzm, 7); + break; + case 7: + tx0 = txm; + ty0 = tym; + tz0 = tzm; + node = node->getChild(7^a); + (st_cur-1)->next = 8; + break; + } + st_cur->next = -1; + } + return nearest_shape; }