pyrit: comparison src/kdtree.cc

equal deleted inserted replaced

-:a5127346fbcd
+:09aedbf5f95f
 			return a.pos < b.pos;
 	};
 };
 // 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):
+	int prev;
-		node(anode), t(at), pb(apb) {};
 };
 // ----------------------------------------
 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;
 	}
 		axis = 2;
 */
 	// 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++)
 		{
 			edges[ax].push_back(ShapeBound(*shape, shapebounds.L[ax], 0));
 	// choose best split pos
 	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++)
 		{
 			if (edge->end)
 	F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl;
 	f += 4;
 #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);
+	getLeftChild()->subdivide(lbb, maxdepth-1);
-	children[1].subdivide(rbb, maxdepth-1);
+	getRightChild()->subdivide(rbb, maxdepth-1);
 }
 void KdTree::build()
 {
 	dbgmsg(1, "* building kd-tree\n");
 	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 */
+	node = root;
-	/* std vector is much faster than stack */
+	StackElem stack[max_depth];
-	vector<StackElem*> st;
+	int entry = 0, exit = 1;
-	StackElem *enPt = new StackElem(NULL, a,
+	stack[entry].t = a;
-		/* distinguish between internal and external origin of a ray*/
-		a >= 0.0 ?
+	/* distinguish between internal and external origin of a ray*/
-			ray.o + ray.dir * a : /* external */
+	if (a >= 0.0)
-			ray.o);               /* internal */
+		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);
+	stack[exit].t = b;
-	st.push_back(exPt);
+	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();
+			splitVal = node->getSplit();
-			short axis = node->getAxis();
+			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;
 				}
 				/* case N4 */
 				farchild = node->getRightChild();
 				node = node->getLeftChild();
 			}
 			else
-			{ /* (enPt->pb[axis] > splitVal) */
+			{ /* (stack[entry].pb[axis] > splitVal) */
-				if (splitVal < exPt->pb[axis])
+				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();
 			}
 			/* 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());
+			register int tmp = exit;
-			exPt->pb.cell[axis] = splitVal;
-			exPt->pb.cell[(axis+1)%3] = ray.o.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
+			exit++;
-			exPt->pb.cell[(axis+2)%3] = ray.o.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
+			if (exit == entry)
-			st.push_back(exPt);
+				exit++;
-		} /* while */
+			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;
-		}
+		entry = exit;
-		enPt = exPt;
+		/* retrieve the pointer to the next node,
-		/* retrieve the pointer to the next node, it is possible that ray traversal terminates */
+		it is possible that ray traversal terminates */
-		node = enPt->node;
+		node = stack[entry].node;
-		st.pop_back();
+		exit = stack[entry].prev;
-	} /* while */
+	}
-	delete enPt;
 	/* ray leaves the scene */
 	return NULL;
 }
 	if (!built)
 		return;
 	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=";
 		save(str, node->getLeftChild());
 		str << "; R=";

branch	pyrit
changeset 74	09aedbf5f95f
parent 72	7c3f38dff082
child 76	3b60fd0bea64