--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/octree.cc Wed Dec 12 19:59:19 2007 +0100
@@ -0,0 +1,325 @@
+/*
+ * 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;
+ BBox shbb;
+ int child, both;
+ unsigned int shapenum = 0;
+ for (sh = shapes->begin(); sh != shapes->end(); sh++)
+ {
+ child = 0;
+ both = 0;
+ shbb = (*sh)->get_bbox();
+
+ if (shbb.L.x >= xsplit)
+ child |= 4; //right
+ else if (shbb.H.x > xsplit)
+ both |= 4; // both
+ // for left, do nothing
+
+ if (shbb.L.y >= ysplit)
+ child |= 2;
+ else if (shbb.H.y > ysplit)
+ both |= 2;
+
+ if (shbb.L.z >= zsplit)
+ child |= 1;
+ else if (shbb.H.z > zsplit)
+ both |= 1;
+
+ if (!both)
+ {
+ getChild(child)->addShape(*sh);
+ shapenum++;
+ }
+ else
+ {
+ // shape goes to more than one child
+ if (both == 7)
+ {
+ for (int i = 0; i < 8; i++)
+ getChild(i)->addShape(*sh);
+ shapenum += 8;
+ }
+ else if (both == 3 || both >= 5)
+ {
+ if (both == 3)
+ {
+ for (int i = 0; i < 4; i++)
+ getChild(child + i)->addShape(*sh);
+ }
+ else if (both == 5)
+ {
+ for (int i = 0; i < 4; i++)
+ getChild(child + i+(i>>1<<1))->addShape(*sh);
+ }
+ else if (both == 6)
+ {
+ for (int i = 0; i < 4; i++)
+ getChild(child + (i<<1))->addShape(*sh);
+ }
+ shapenum += 4;
+ }
+ else
+ {
+ getChild(child)->addShape(*sh);
+ getChild(child+both)->addShape(*sh);
+ shapenum += 2;
+ }
+ }
+ }
+
+ 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;
+
+ //Octree *node = root;
+ 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;
+}