--- a/src/kdtree.cc	Sat Apr 19 18:00:27 2008 +0200
+++ b/src/kdtree.cc	Sun Apr 20 16:48:24 2008 +0200
@@ -78,60 +78,65 @@
 	}
 
 	// choose split axis
-	axis = 0;
+	/*axis = 0;
 	if (bounds.h() > bounds.w() && bounds.h() > bounds.d())
 		axis = 1;
 	if (bounds.d() > bounds.w() && bounds.d() > bounds.h())
 		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++)
 	{
 		BBox shapebounds = (*shape)->get_bbox();
-	//	for (int ax = 0; ax < 3; ax++)
-	//	{
-			edges[axis].push_back(ShapeBound(*shape, shapebounds.L[axis], 0));
-			edges[axis].push_back(ShapeBound(*shape, shapebounds.H[axis], 1));
-	//	}
+		for (int ax = 0; ax < 3; ax++)
+		{
+			edges[ax].push_back(ShapeBound(*shape, shapebounds.L[ax], 0));
+			edges[ax].push_back(ShapeBound(*shape, shapebounds.H[ax], 1));
+		}
 	}
-	sort(edges[axis].begin(), edges[axis].end());
+	for (int ax = 0; ax < 3; ax++)
+		sort(edges[ax].begin(), edges[ax].end());
 
 	// 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
-	BBox lbb = bounds;
-	BBox rbb = bounds;
 
-	vector<ShapeBound>::iterator edge, splitedge = edges[axis].end();
-	int lnum = 0, rnum = shapes->size();
-	for (edge = edges[axis].begin(); edge != edges[axis].end(); edge++)
+	vector<ShapeBound>::iterator edge, splitedge = edges[2].end();
+	for (int ax = 0; ax < 3; ax++)
 	{
-		if (edge->end)
-			rnum--;
+		int lnum = 0, rnum = shapes->size();
+		BBox lbb = bounds;
+		BBox rbb = bounds;
+		for (edge = edges[ax].begin(); edge != edges[ax].end(); edge++)
+		{
+			if (edge->end)
+				rnum--;
 
-		// calculate SAH cost of this split
-		lbb.H.cell[axis] = edge->pos;
-		rbb.L.cell[axis] = edge->pos;
-		Float SAL = (lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
-		Float SAR = (rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
-		Float splitcost = K*SAV + SAL*(K + lnum) + SAR*(K + rnum);
+			// calculate SAH cost of this split
+			lbb.H.cell[ax] = edge->pos;
+			rbb.L.cell[ax] = edge->pos;
+			Float SAL = (lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
+			Float SAR = (rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
+			Float splitcost = K*SAV + SAL*(K + lnum) + SAR*(K + rnum);
 
-		if (splitcost < cost)
-		{
-			splitedge = edge;
-			cost = splitcost;
-			split = edge->pos;
+			if (splitcost < cost)
+			{
+				axis = ax;
+				splitedge = edge;
+				cost = splitcost;
+				split = edge->pos;
+			}
+
+			if (!edge->end)
+				lnum++;
 		}
-
-		if (!edge->end)
-			lnum++;
 	}
 
-	if (splitedge == edges[axis].end())
+	if (splitedge == edges[2].end())
 	{
 		setLeaf();
 		return;
@@ -162,17 +167,19 @@
 
 	// split this node
 	delete shapes;
+	BBox lbb = bounds;
+	BBox rbb = bounds;
+	lbb.H.cell[axis] = split;
+	rbb.L.cell[axis] = split;
 	children = new KdNode[2];
+
 	for (edge = edges[axis].begin(); edge != splitedge; edge++)
-		if (!edge->end)
+		if (!edge->end && edge->shape->intersect_bbox(lbb))
 			children[0].addShape(edge->shape);
 	for (edge = splitedge; edge < edges[axis].end(); edge++)
-		if (edge->end)
+		if (edge->end && edge->shape->intersect_bbox(rbb))
 			children[1].addShape(edge->shape);
 
-	lbb.H.cell[axis] = split;
-	rbb.L.cell[axis] = split;
-
 	children[0].subdivide(lbb, maxdepth-1);
 	children[1].subdivide(rbb, maxdepth-1);
 }