#if TOIMPLEMENT

#include <queue>
#include <stack>
#include "Bvh.h"
#include "yare.h"
#include "PerfTimer.h"
#include "Camera.h"
#include "Settings.h"
#include "Context.h"
#include "TriangleBvh.h"
#include "NodeGeometry.h"
#include "Viewer.h"


#include <fstream>
#include <iostream>
#include <iomanip>

#define INVALID_TEST ((unsigned int)-1)
#define MAX_FLOAT 1e20f

#define TYPE_INTERIOR -2
#define TYPE_LEAF -3
#define USE_TIGHTER_BOUNDS_FOR_ALL 0

const static bool useVbos = true;
unsigned int Bvh::sCurrentVboId = -1;

using namespace std;



/*
3 x---------x 2
  |\         \
  | \         \
  |  \         \
  |	4 x---------x 5
  |   |         | 
0 x   |     x 1 |
   \  |         |
    \ |         |
     \|         |
  	7 x---------x 6             
*/

static unsigned int sIndices[] =	
{7, // degenerated
 7, 6, 4, 5, 3, 2, 0, 1,
 1, 4, // degenerated
 4, 3, 7, 0, 6, 1, 5, 2,
 2 // degenerated
};


const static int sNumIndicesPerBox = 20;

/*	Order of vertices
	0 = (0, 0, 0)
	1 = (1, 0, 0)
	2 = (1, 1, 0)
	3 = (0, 1, 0)
	4 = (0, 1, 1)
	5 = (1, 1, 1)
	6 = (1, 0, 1)
	7 = (0, 0, 1)
*/


BvhNode::BvhNode(BvhNode *parent): 
mParent(parent),
mAxis(-1), 
mDepth(0), 
mPlaneMask(0),
mPreferredPlane(0),
mVizBox(NULL),
mLastRenderedFrame(-999),
mFirst(-1), 
mLast(-1),
mNumTestNodes(1),
mTestNodesIdx(-1),
mIndicesPtr(-1)
{
}


void BvhNode::ResetVisibility()
{
	mVisibility.Reset();
	mMeasurements.Reset();

	mLastRenderedFrame = -999;
}


BvhNode::~BvhNode() 
{
	if (mVizBox) mVizBox->unref();
}


Box *BvhNode::GetOrCreateVizBox() 
{
	if (!mVizBox) 
	{
		mVizBox = new Box(mBox); 
		mVizBox->ref();
	}

	return mVizBox;
}


void BvhNode::VisibilityInfo::Reset()
{
	mIsVisible = false;
	mIsFullyVisible = false;

	mVisiblePixels = 0;
	mAssumedVisibleFrames = 0;
	mRemainingVisibleFrames = 0;

	mLastVisitedFrame = -1;
	mLastTestedFrame = -1;
	mTurnedVisibleFrame = 0;

	mTimesInvisible = 0;
	mTimesTested = 0;
	mTimesChangedClassification = 0;
	mAvgChangedClassification = 0;
	mIsFrustumCulled = false;

	mLastTestedVisibleFrame = 0;
	mIsNew = true;
}


void BvhNode::SetFullyVisible(const bool visible)
{
	mVisibility.mIsFullyVisible = visible;
}


BvhLeaf::~BvhLeaf()
{
	if (mTriangleBvh) delete mTriangleBvh;
}


/***********************************************************/
/*                class Bvh implementation                 */
/***********************************************************/



Bvh::Bvh(const GeometryVector &geometry, 
		 DistanceSortedRenderAction *const renderer):
mCamera(NULL), 
mFrameId(-1), 
mRoot(NULL), 
mVertices(NULL),
mIndices(NULL),
mUseTighterBoundsOnlyForLeafTests(false),
mRenderer(renderer),
mTestIndices(NULL),
mCurrentIndicesPtr(0),
mCollectTighterBoundsWithMaxLevel(true)
{
	mGeometry = new NodeGeometry*[geometry.size()];
	mGeometrySize = geometry.size();

	BoundingBox sceneBox;

	// compute scene extent
	for (size_t i = 0; i < mGeometrySize; ++ i)
	{
		mGeometry[i] = geometry[i];
		sceneBox.combine(&mGeometry[i]->mBox);

		mBvhStats.mTriangles += mGeometry[i]->mNumTriangles;
	}


	////////
	//-- create root

	BvhLeaf *leaf = new BvhLeaf(NULL);
	leaf->mDepth = 0;
	leaf->mFirst = 0;
	leaf->mLast = (int)geometry.size() - 1;
	leaf->mBox = sceneBox;

	OUT1("geometry in root: " << leaf->mLast);

	mRoot = leaf;
	mNumNodes = 1;

	// parameters
	mMaxGeometry = Settings::Global()->get_nvocc_bvh_max_objects();
	mMaxTriangles = Settings::Global()->get_nvocc_bvh_max_triangles();
	mMaxDepth = Settings::Global()->get_nvocc_bvh_max_depth();
	mSplitType = Settings::Global()->get_nvocc_bvh_split_type();

	float sceneArea = sceneBox.getSurface();
	float minAreaRatio = Settings::Global()->get_nvocc_bvh_min_area();

	mMinArea = minAreaRatio * sceneArea;

	mMaxDepthForTestingChildren = Settings::Global()->get_nvocc_bvh_max_depth_for_testing_children();
	mUseTighterBoundsOnlyForLeafTests = (Settings::Global()->get_nvocc_use_tighter_bounds() >= 1);
	
	mAreaRatioThreshold = Settings::Global()->get_nvocc_area_ratio_threshold();
	mVolRatioThreshold = Settings::Global()->get_nvocc_vol_ratio_threshold();
}


Bvh::~Bvh() 
{
	if (mVertices) delete []mVertices; 
	if (mIndices) delete [] mIndices;
	if (mTestIndices) delete [] mTestIndices;
	if (mGeometry) delete [] mGeometry;

	if (mRoot) delete mRoot;
}


void Bvh::Construct()
{
	PerfTimer constructTimer;
	constructTimer.Entry();

	OUT1("Info: Constructing BVH");

	mRoot = SubdivideLeaf((BvhLeaf *)mRoot, 0);

	constructTimer.Exit();

	OUT1("Info: Bvh done in " << constructTimer.TotalCount() << " seconds.");
	OUT1("Info: Number Of BVH Nodes = " << mNumNodes);

	// update stats on leaf nodes
	UpdateNumLeaves(mRoot);
	
	mAvgDepth = 1.0f + log((float)GetNumLeaves()) / log(2.0f);

	BvhLeafContainer leaves;
	leaves.reserve(GetNumLeaves());
	CollectLeaves(mRoot, leaves);

	// compute tighter boundaries for the leaves
	PostProcessLeaves(leaves);
	// compute unique ids
	ComputeIds();
	/// pull up some data from the leaves
	UpdateInteriors(mRoot);
	// recompute the boundaries once
	RecomputeBounds();
	// do stats
	ComputeBvhStats();

	PrintBvhStats();
}


float Bvh::EvaluateSahCost(BvhLeaf *leaf, const int axis, const float position)
{
	// count triangles on the left / right
	int left = 0;
	int right = 0;
	BoundingBox leftBox, rightBox;

	const int candidates = std::max(50, (int)(leaf->CountPrimitives() / 20));

	const float finc = leaf->CountPrimitives() / (float)candidates;

	int i = leaf->mFirst;
	float fi = leaf->mFirst + 0.5f;

	BoundingBox box;

	for (; i <= leaf->mLast;) 
	{
		box = mGeometry[i]->GetBoundingBox();

		if (box.getCenter()[axis] < position) 
		{
			++ left;
			// update the box
			leftBox.combine(&box);
		} 
		else 
		{
			++ right;
			rightBox.combine(&box);
		}
		 
		fi += finc;
		i = (int)fi;
	}

	float bW = 1.0f;
	float leftRatio = left / (float)leaf->CountPrimitives();
	float rightRatio = right / (float)leaf->CountPrimitives();

	float saLeft = 0.0f;
	float saRight = 0.0f;

	// not a valid split
	if (!left || !right)
		return 1e25f;

	saLeft = leftBox.getSurface();
	saRight = rightBox.getSurface();


	return
		saLeft  * ((1.0f - bW) + bW * leftRatio) + 
		saRight * ((1.0f - bW) + bW * rightRatio);
}


float Bvh::SelectPlaneSah(BvhLeaf *leaf, int &axis, float &minCost)
{
	minCost = MAX_FLOAT;
	float bestPos = minCost;
	int bestAxis = 0;

	//  cout<<"Evaluating axis..."<<endl<<flush;

	const int initialPlanes = 3;

	// initiate the costs
	for (axis = 0; axis < 3; ++ axis) 
	{
		const float size = leaf->mBox.getUpper()[axis] - leaf->mBox.getLower()[axis];

		for (int i = 0; i < initialPlanes; ++ i) 
		{
			const float pos = leaf->mBox.getLower()[axis] + (i + 1) * size / (initialPlanes + 1);
			const float cost = EvaluateSahCost(leaf, axis, pos);

			if (cost < minCost) 
			{
				minCost = cost;
				bestPos = pos;
				bestAxis = axis;
			}
		}
	}

	axis = bestAxis;

	//  cout<<axis<<endl<<flush;
	const float shrink = 0.5f;

	// now hierarchically refine the initial interval
	float size = shrink * 
		(leaf->mBox.getUpper()[axis] - leaf->mBox.getLower()[axis]) / (initialPlanes + 1);

	const int steps = 4;
	float cost;

	for (int i = 0; i < steps; ++ i) 
	{
		const float left = bestPos - size;
		const float right = bestPos + size;

		cost = EvaluateSahCost(leaf, axis, left);

		if (cost < minCost) 
		{
			minCost = cost;
			bestPos = left;
		}

		cost = EvaluateSahCost(leaf, axis, right);

		if (cost < minCost) 
		{
			minCost = cost;
			bestPos = right;
		}

		size = shrink * size;
	}

	//OUT1("best axis: " << axis << " " << bestPos << " " << minCost);

	return bestPos;
}


void Bvh::ResizeVisibilityBuffers(const int maxSize)
{
	std::stack<BvhNode *> nodeStack;

	nodeStack.push(mRoot);

	while (!nodeStack.empty()) 
	{
		BvhNode *node = nodeStack.top();
		nodeStack.pop();

		static_cast<BvhLeaf *>(node)->mMeasurements.Reset();

		if (!node->IsLeaf()) 
		{
			BvhInterior *interior = static_cast<BvhInterior *>(node);
			
			nodeStack.push(interior->mBack);
			nodeStack.push(interior->mFront);
		}
	}
}


int Bvh::SortTriangles(BvhLeaf *leaf, 
					   const int axis, 
					   const float position)
{
	int i = leaf->mFirst;
	int j = leaf->mLast;

    while (1) 
	{
		while (mGeometry[i]->GetBoundingBox().getCenter()[axis] < position) 
			++ i;
	
		while (position < mGeometry[j]->GetBoundingBox().getCenter()[axis])
			-- j;

		if (i < j) 
		{
			std::swap(mGeometry[i], mGeometry[j]);
			
			++ i;
			-- j;
		}
		else
		{
			break;
		}
	}

	return j;
}


int Bvh::SortTrianglesSpatialMedian(BvhLeaf *leaf, const int axis)
{
	// spatial median
	float x = leaf->mBox.getCenter()[axis];

	return SortTriangles(leaf, axis, x);
}


int Bvh::SortTrianglesObjectMedian(BvhLeaf *leaf, const int axis, float &pos)
{
	// Now distribute the objects into the subnodes
	// Use QuickMedian sort
	int l = leaf->mFirst;
	int r = leaf->mLast;
	int k = (l + r) / 2;

	while (l < r) 
	{
		int i = l;
		int j = r;

		// get some estimation of the pivot
		pos = mGeometry[(l + r) / 2]->GetBoundingBox().getCenter()[axis];

		while (1) 
		{
			while ((i  <= leaf->mLast) && (mGeometry[i]->GetBoundingBox().getCenter()[axis] < pos))
				++ i;

			while((j >= leaf->mFirst) && (pos < mGeometry[j]->GetBoundingBox().getCenter()[axis]))
				-- j;

			if (i <= j) 
			{
				std::swap(mGeometry[i], mGeometry[j]);
				++ i;
				-- j;
			}
			else
			{
				break;
			}
		}

		// now check the extents
		if (i >= k)
			r = j;
		else
			l = i;
	}

	return k;
}


int Bvh::SortTrianglesSurfaceArea(BvhLeaf *leaf, const float sa)
{
	int i = leaf->mFirst;
	int j = leaf->mLast;

	while(1) 
	{
		while ((i <= j) && (mGeometry[i]->GetBoundingBox().getSurface() < sa)) 
			++ i;

		while ((i <= j) && (sa < mGeometry[j]->GetBoundingBox().getSurface()))
			-- j;

		if (i < j) 
		{
			swap(mGeometry[i], mGeometry[j]);
			++ i;
			-- j;
		}
		else
			break;
	}

	return j;
}


bool Bvh::TerminationCriteriaMet(BvhLeaf *leaf) const
{
	const bool terminationCriteriaMet = 
		(leaf->mBox.getSurface() < mMinArea) ||
		(leaf->CountPrimitives() <= mMaxGeometry) ||
		(CountTriangles(leaf) <= mMaxTriangles) ||
		(leaf->mDepth > mMaxDepth);

	return terminationCriteriaMet;
}


BvhNode *Bvh::SubdivideLeaf(BvhLeaf *leaf, const int parentAxis)
{
	if (TerminationCriteriaMet(leaf))
		return leaf;
	
	//int axis = leaf->mBox.MajorAxis();
	int axis = (parentAxis + 1) % 3;

	// position of the split in the partailly sorted array of triangles
	// corresponding to this leaf
	int split = -1;
	const int scale = 20;

	float pos = -1.0f;

	// Spatial median subdivision
	switch (mSplitType) 
	{
	case SPATIAL_MEDIAN:
		split = SortTrianglesSpatialMedian(leaf, axis);

		if ( 
			((split - leaf->mFirst) < leaf->CountPrimitives() / scale) ||
			((leaf->mLast - split) < leaf->CountPrimitives() / scale) ) 
		{
			split = SortTrianglesObjectMedian(leaf, axis, pos);
		}

		pos = leaf->mBox.getCenter()[axis];
		break;
	case OBJECT_MEDIAN:
		// Object median subdivision: approximately the same number
		// of objects on the left of the the splitting point.
		split = SortTrianglesObjectMedian(leaf, axis, pos);
		break;
	case SAH:
		{
			float cost;
			pos = SelectPlaneSah(leaf, axis, cost);

			if (pos != MAX_FLOAT)
			{
				split = SortTriangles(leaf, axis, pos);
			}

			if ((pos == MAX_FLOAT) || (split == leaf->mLast))
			{
				split = -1;
				split = SortTrianglesObjectMedian(leaf, axis, pos);
			}
		}
		break;
	case SAH_OR_SIZE: 
		{
			// split by size instead
			const float saThreshold = 0.2f * leaf->GetBox().getSurface();
			split = SortTrianglesSurfaceArea(leaf, saThreshold);

			if ((split == leaf->mLast) || (split == leaf->mFirst - 1)) 
			{
				// use SAH
				float cost;
				pos = SelectPlaneSah(leaf, axis, cost);

				if (pos != MAX_FLOAT)
					split = SortTriangles(leaf, axis, pos);
				else
					split = SortTrianglesObjectMedian(leaf, axis, pos);
			}
			else
			{
				// note: no position is computed!!
				//OUT1("sorted by size");
			}
		}
		break;
	default:
		OUT1("should not come here");
		break;
	}

	if (1 && ((split == leaf->mLast)))
	{
		// no reduction: we should never come here
		OUT1("error: no reduction " << leaf->CountPrimitives() << " " << leaf->mFirst << " " << split << " " << leaf->mLast);
		return leaf;
	}

	// create two more nodes
	mNumNodes += 2;
	BvhInterior *parent = new BvhInterior(leaf->GetParent());
	BvhLeaf *front = new BvhLeaf(parent);
	
	parent->mAxis = axis;
	parent->mBox = leaf->mBox;
	parent->mDepth = leaf->mDepth;

	parent->mBack = leaf;
	parent->mFront = front;
	parent->mPosition = pos;
	
	// now assign the triangles to the subnodes
	front->mFirst = split + 1;
	front->mLast = leaf->mLast;
	front->mDepth = leaf->mDepth + 1;
	leaf->mLast = split;
	leaf->mDepth = front->mDepth;
	leaf->mParent = parent;
	
	// reset box
	leaf->mBox = BoundingBox();

	UpdateLeafBox(static_cast<BvhLeaf *>(parent->mBack));
	UpdateLeafBox(static_cast<BvhLeaf *>(parent->mFront));

	// some output
	const int n = 500;
	if ((GetNumLeaves() % n) == n - 1)
		OUT1("created " << GetNumLeaves() << " leaves");
	
#if VERBOSE_OUTPUT
	OUT1("box: " << parent->mBox.getUpper() << " " << parent->mBox.getLower());
	OUT1("depth: " << (int)parent->mDepth);
	OUT1("axis: " << axis);
	OUT1("bc="<<((BvhLeaf *)parent->mBack)->Count());
	OUT1("fc="<<((BvhLeaf *)parent->mFront)->Count());
	OUT1("back: " << parent->mBack->mBox.getSurface());
	OUT1("front: " << parent->mFront->mBox.getSurface());
#endif

	// recursively continue subdivision
	parent->mBack = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mBack), axis);
	parent->mFront = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mFront), axis);

	return parent;
}


void Bvh::UpdateLeafBox(BvhLeaf *leaf)
{
	for (int i = leaf->mFirst; i <= leaf->mLast; ++ i) 
	{
		leaf->mBox.combine(&mGeometry[i]->mBox);
	} // for
}


void Bvh::UpdateNumLeaves(BvhNode *node) const
{
	if (node->IsLeaf())
	{
		node->mNumLeaves = 1;
	}
	else
	{
		BvhNode *f = static_cast<BvhInterior *>(node)->GetFront();
		BvhNode *b = static_cast<BvhInterior *>(node)->GetBack();

		UpdateNumLeaves(f);
		UpdateNumLeaves(b);

		node->mNumLeaves = f->mNumLeaves + b->mNumLeaves;
	}
}


void Bvh::UpdateBoxes(BvhNode *node)
{
	if (node->IsLeaf()) 
	{
		UpdateLeafBox(static_cast<BvhLeaf *>(node));
	}
	else 
	{
		BvhInterior *interior = static_cast<BvhInterior *>(node);
	
		UpdateBoxes(interior->mBack);
		UpdateBoxes(interior->mFront);

		node->mBox = interior->mBack->mBox;
		node->mBox.combine(&interior->mFront->mBox);
	}
}


void Bvh::UpdateFullVisibility(BvhNode *node) const
{		
	// node not visited in this frame => no change
	if (node->GetLastVisitedFrame() != mFrameId) 
		return;

	// leaf node: terminate recursion
	if (node->IsLeaf())
	{
		node->SetFullyVisible(node->IsVisible());	
		return;
	}

	BvhInterior *interior = static_cast<BvhInterior *>(node);
	
	// recursive traversal
	UpdateFullVisibility(interior->mBack);
	UpdateFullVisibility(interior->mFront);
	
	interior->SetFullyVisible(interior->mBack->IsFullyVisible() && 
							  interior->mFront->IsFullyVisible());
}


void Bvh::PullUpLastVisited(BvhNode *node, const int frameId) const
{		
	BvhNode *parent = node->GetParent();

	while (parent && (parent->GetLastVisitedFrame() != frameId))
	{
		parent->SetLastVisitedFrame(frameId);
		parent = parent->GetParent();
	}
}


void Bvh::MakeParentsVisible(BvhNode *node)
{
	BvhNode *parent = node->GetParent();

	while (parent && (!parent->IsVisible()))
	{
		parent->SetVisible(true);
		parent = parent->GetParent();
	}
}


void Bvh::CollectLeaves(BvhNode *node, BvhLeafContainer &leaves)
{
	stack<BvhNode *> tStack;
	tStack.push(node);

	while (!tStack.empty())
	{
		BvhNode *node = tStack.top();

		tStack.pop();

		if (!node->IsLeaf())
		{
			BvhInterior *interior = static_cast<BvhInterior *>(node);

			tStack.push(interior->mFront);
			tStack.push(interior->mBack);
		}
		else
		{
			leaves.push_back(static_cast<BvhLeaf *>(node));
		}
	}
}


void Bvh::CollectNodes(BvhNode *node, BvhNodeContainer &nodes)
{
	stack<BvhNode *> tStack;

	tStack.push(node);

	while (!tStack.empty())
	{
		BvhNode *node = tStack.top();
		tStack.pop();

		nodes.push_back(node);

		if (!node->IsLeaf())
		{
			BvhInterior *interior = static_cast<BvhInterior *>(node);

			tStack.push(interior->mFront);
			tStack.push(interior->mBack);
		}
	}
}


typedef pair<BvhNode *, int> tPair;

void Bvh::CollectNodes(BvhNode *root, 
					   const int depth, 
					   HierarchyNodeContainer &nodes)
{
	stack<tPair> tStack;
	tStack.push(tPair(root, 0));

	while (!tStack.empty())
	{
		BvhNode *node = tStack.top().first;
		const int d = tStack.top().second;

		tStack.pop();

		// found depth => take this node
		if (d == depth)
		{
			nodes.push_back(node);
		}
		else if (node->IsLeaf())
		{
			// found leaf
			BvhLeaf *leaf = static_cast<BvhLeaf *>(node);

			// there is no further subdivision on triangle level
			if (leaf->mTriangleBvh->GetNumNodes() == 1)
			{
				nodes.push_back(leaf);
			}
			else // more than a root => search in local bvh
			{
				leaf->mTriangleBvh->
					CollectNodes(leaf->mTriangleBvh->GetRoot(), depth - d, nodes);
			}
		}
		else // interior
		{
			BvhInterior *interior = static_cast<BvhInterior *>(node);

			tStack.push(tPair(interior->mFront, d + 1));
			tStack.push(tPair(interior->mBack, d + 1));
		}
	}
}


void Bvh::CollectAllNodes(BvhNode *node, HierarchyNodeContainer &nodes)
{
	BvhNodeContainer bvhNodes;

	// collect nodes of geometry bvh
	CollectNodes(mRoot, bvhNodes);

	
	// first collect nodes of object bvh
	BvhNodeContainer::const_iterator lit, lit_end = bvhNodes.end();

	for (lit = bvhNodes.begin(); lit != lit_end; ++ lit)
	{
		nodes.push_back(*lit);
	}

	// collect nodes of local bvh: all except roots, 
	// because they are equivalent to geometry bvh leaves
	for (lit = bvhNodes.begin(); lit != lit_end; ++ lit)
	{
		if ((*lit)->IsLeaf())
		{
			BvhLeaf *leaf = static_cast<BvhLeaf *>(*lit);

			TriangleBvhNodeContainer hnodes;
			leaf->mTriangleBvh->CollectNodes(leaf->mTriangleBvh->GetRoot(), hnodes);

			TriangleBvhNodeContainer::const_iterator hit, hit_end = hnodes.end();

			for (hit = hnodes.begin(); hit != hit_end; ++ hit)
			{
				TriangleBvhNode *n = *hit;

				// don't include root of local bvh - it has the same bounds as the leaves
				if (n != leaf->mTriangleBvh->GetRoot())
					nodes.push_back(n);
			}
		}
	}
}


void Bvh::CollectVisibleLeaves(BvhNode *node, BvhLeafContainer &leaves)
{
	stack<BvhNode *> tStack;
	tStack.push(node);

	while (!tStack.empty())
	{
		BvhNode *node = tStack.top();
		tStack.pop();

		if (IsWithinViewFrustum(node))
		{
			if (!node->IsLeaf())
			{
				BvhInterior *interior = static_cast<BvhInterior *>(node);

				tStack.push(interior->mFront);
				tStack.push(interior->mBack);
			}
			else
			{
				leaves.push_back(static_cast<BvhLeaf *>(node));
			}
		}
	}
}


BvhLeaf *Bvh::GetRandomLeaf(BvhNode *node)
{
	stack<BvhNode *> nodeStack;
  
	nodeStack.push(node);

	int mask = rand();

    while (!nodeStack.empty()) 
	{
		BvhNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf()) 
		{
			// random leaf
			return static_cast<BvhLeaf *>(node);
		} 
		else 
		{
			BvhInterior *interior = static_cast<BvhInterior *>(node);
			
			BvhNode *next;
			
			// random decision
			if (mask & 1)
				next = interior->mBack;
			else
				next = interior->mFront;

			mask = mask >> 1;

			nodeStack.push(next);
		}
	}
  
	// should never come here
	return NULL;
}


BvhLeaf *Bvh::GetRandomVisibleLeaf(BvhNode *node)
{
	BvhLeafContainer leaves;
	leaves.reserve(node->GetNumLeaves());

	CollectVisibleLeaves(node, leaves);

	if (leaves.empty())
		return NULL;

	const int r = (int)(((float)rand() / RAND_MAX) * ((float)leaves.size() - 0.5f));

	return leaves[r];
}


void Bvh::CollectGeometry(BvhNode *node, GeometryVector &geometry)
{
	//geometry.reserve(node->CountPrimitives());

	for (int i = node->mFirst; i <= node->mLast; ++ i)
	{
		geometry.push_back(mGeometry[i]);
	}
}


NodeGeometry **Bvh::GetGeometry(BvhNode *node, int &geometrySize)
{
	geometrySize = node->CountPrimitives();
	return mGeometry + node->mFirst;
}


void Bvh::UpdateInteriors(BvhNode *node)
{
	if (!node->IsLeaf())
	{
		BvhInterior *interior = static_cast<BvhInterior *> (node);
		
		// update the indices of the geometry so we can render interiors as well
		UpdateInteriors(interior->GetBack());
		UpdateInteriors(interior->GetFront());

		// update area
		interior->mFirst = min(interior->GetBack()->mFirst, interior->GetFront()->mFirst);
		interior->mLast = max(interior->GetBack()->mLast, interior->GetFront()->mLast);
	
		interior->mArea = interior->GetBack()->mArea + interior->GetFront()->mArea;
	}
}


int	Bvh::IsWithinViewFrustumLocal(BvhNode *node)
{
	bool bIntersect	= false;

	if (node->GetParent())
		node->mPlaneMask = node->GetParent()->mPlaneMask;


	////////
	//-- do the view frustum culling for the planes [mPreferredPlane - 5]

	for (int i = node->mPreferredPlane; i < 6; ++ i)
	{
		//-- do the test only if necessary
		if (node->mPlaneMask & (1 << i))
		{
			//-- test the n-vertex
			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
			{
				//-- outside
				node->mPreferredPlane = i;
				return 0;
			}

			//-- test the p-vertex
			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f)
			{
				//-- completely inside: children don't need to check against this plane no more
				node->mPlaneMask^=	1 << i;
			}
			else 
			{
				bIntersect = true;
			}
		}
	}

	//////////
	//-- do the view frustum culling for the planes [0 - m_iPreferredPlane)

	for (int i = 0; i < node->mPreferredPlane; ++ i)
	{
		// do the test only if necessary
		if (node->mPlaneMask & (1 << i))
		{
			//-- test the n-vertex
			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
			{
				// outside
				node->mPreferredPlane = i;
				return 0;
			}

			//-- test the p-vertex
			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f)
			{
				// completely inside: children don't need to check against this plane no more
				node->mPlaneMask^= 1 << i;
			}
			else 
			{
				bIntersect = true;
			}
		}
	}

	return bIntersect ? -1 : 1;
}


int	Bvh::IsWithinViewFrustum(const BoundingBox &box, const int planeMask, const int preferredPlane)
{
	bool bIntersect	= false;

	////////
	//-- do the view frustum culling for the planes [mPreferredPlane - 5]

	for (int i = preferredPlane; i < 6; ++ i)
	{
		//-- do the test only if necessary
		if (planeMask & (1 << i))
		{
			//-- test the n-vertex
			if (box.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
				return 0;

			//-- test the p-vertex
			if (!(box.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f))
				bIntersect = true;
		}
	}

	//////////
	//-- do the view frustum culling for the planes [0 - m_iPreferredPlane)

	for (int i = 0; i < preferredPlane; ++ i)
	{
		// do the test only if necessary
		if (planeMask & (1 << i))
		{
			//-- test the n-vertex
			if (box.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
				// outside
				return 0;

			//-- test the p-vertex
			if (!(box.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f))
				bIntersect = true;
		}
	}

	return bIntersect ? -1 : 1;
}


int	Bvh::IsWithinViewFrustum(BvhNode *node)
{
	if (node->mCache.mLastFrustumTestedFrameId == mFrameId)
	{
		return node->mCache.mFrustumCulled;
	}

	node->mCache.mLastFrustumTestedFrameId = mFrameId;

	timeViewFrustumCulling.Entry();
	
	if (Settings::Global()->get_nvocc_use_tighter_bounds_for_frustum_culling())
	{
		const int intersect = IsWithinViewFrustumLocal(node);

		if ((node->mNumTestNodes == 1) || (intersect != -1))
		{
			//OUT1("x");
			node->mCache.mFrustumCulled = intersect;
		}
		else 
		{
			// maybe intersecting one of the tighter boxes
			node->mCache.mFrustumCulled = 0;

			for (int i = 0; i < node->mNumTestNodes; ++ i)
			{
				RenderableHierarchyNode *n = mTestNodes[node->mTestNodesIdx + i];

				if (IsWithinViewFrustum(n->GetBox(), 
										node->mPlaneMask, 
										node->mPreferredPlane))
				{
					node->mCache.mFrustumCulled = -1;
					break;
				}
			}
		}
	}
	else
	{
		//OUT1("y");
		node->mCache.mFrustumCulled = IsWithinViewFrustumLocal(node);
	}

	timeViewFrustumCulling.Exit();

	return node->mCache.mFrustumCulled;
}


void Bvh::InitFrame(Camera *camera, const int currentFrameId, Viewer *viewer)
{
	// = 0011 1111 which means that at the beginning, all six planes have to frustum culled
	mRoot->mPlaneMask = 0x3f;
	mCamera = camera;

	mFrameId = currentFrameId;

	// to begin with, we must grab the plane equations of the six clipplanes of the viewfrustum
	Matrix4f matViewing, matProjectionView;
	Transform3D	transform;

	camera->GetTransform(transform);
	transform.get(matViewing);

	camera->GetProjectionMatrix(matProjectionView);
	matProjectionView *= matViewing;

	Vec3f vec;
	float fInvLength;


	//////////
	//-- extract the plane equations

	for (int i = 0; i < 4; ++ i)
	{
		mFrustum.mClipPlane[0][i] = matProjectionView[i][3] - matProjectionView[i][0];	// right plane
		mFrustum.mClipPlane[1][i] = matProjectionView[i][3] + matProjectionView[i][0];	// left plane
		mFrustum.mClipPlane[2][i] = matProjectionView[i][3] + matProjectionView[i][1];	// bottom plane
		mFrustum.mClipPlane[3][i] = matProjectionView[i][3] - matProjectionView[i][1];	// top plane
		mFrustum.mClipPlane[4][i] = matProjectionView[i][3] - matProjectionView[i][2];	// far plane
		mFrustum.mClipPlane[5][i] = matProjectionView[i][3] + matProjectionView[i][2];	// near plane
	}

	////////////
	//-- normalize the coefficients and find the indices of the n- and p-vertices

	for (int i = 0; i < 6; ++ i)
	{
		// the clipping planes look outward the frustum,
		// so distances > 0 mean that a point is outside
		fInvLength	= -1.0f / sqrt(	mFrustum.mClipPlane[i][0] * mFrustum.mClipPlane[i][0] +	
									mFrustum.mClipPlane[i][1] * mFrustum.mClipPlane[i][1] +	
									mFrustum.mClipPlane[i][2] * mFrustum.mClipPlane[i][2]);

		mFrustum.mClipPlane[i][0] *= fInvLength;
		mFrustum.mClipPlane[i][1] *= fInvLength;
		mFrustum.mClipPlane[i][2] *= fInvLength;
		mFrustum.mClipPlane[i][3] *= fInvLength;

		vec.x = mFrustum.mClipPlane[i][0];
		vec.y = mFrustum.mClipPlane[i][1];
		vec.z = mFrustum.mClipPlane[i][2];

		// n-vertex
		mClipPlaneAABBVertexIndices[i][0] = BoundingBox::getIndexNearestVertex(vec);
		// p-vertex
		mClipPlaneAABBVertexIndices[i][1] = BoundingBox::getIndexFurthestVertex(vec);	
	}


	const float aspect = mCamera->GetAspect();
	const float fov = mCamera->GetFov();
	
	mScale = sqrt(aspect) * 2.0f * tan(fov * mypi / 180.0f);

	mNearPlane = mCamera->GetRealViewDirection();
	mNearPlaneD = - mNearPlane * (mCamera->GetPosition() - Origin3f);
}


float Bvh::CalcDistance(BvhNode *node) const
{
	timeDistance.Entry();

	float distance;

#if USE_TIGHTER_BOUNDS_FOR_ALL
	float minDist = 1e25f;

	//OUT1("testing " << mNumTestNodes << " nodes");
	for (int i = 0; i < node->mNumTestNodes; ++ i)
	{
		RenderableHierarchyNode *hnode = mTestNodes[node->mTestNodesIdx + i];

		if(hnode->mCache.mLastDistanceTestedFrameId < mFrameId)
		{
			hnode->mCache.mDistance = 
				hnode->GetBox().getMinVisibleDistance(mNearPlane, mNearPlaneD);
			hnode->mCache.mLastDistanceTestedFrameId = mFrameId;
		}

		if (hnode->mCache.mDistance < minDist)
			minDist = hnode->mCache.mDistance;
	}

	distance = minDist;

#else

	if(node->mCache.mLastDistanceTestedFrameId != mFrameId)
	{
		//OUT1("z " << node->mLastDistanceTestedFrameId  << " " << mFrameId);
		node->mCache.mDistance = node->GetBox().getMinVisibleDistance(mNearPlane, mNearPlaneD);
		node->mCache.mLastDistanceTestedFrameId = mFrameId;
	}

	distance = node->mCache.mDistance;
#endif

	timeDistance.Exit();

	return distance;
}


float Bvh::GetSquareDistance(BvhNode *node) const
{
	timeDistance.Entry();

	float distance;

#if USE_TIGHTER_BOUNDS_FOR_ALL
	const Vector3f nearplane = mCamera->GetRealViewDirection();
	const float nearplaneD = - nearplane * (mCamera->GetPosition() - Origin3f);

	float minDist = 1e25f;
	
	//OUT1("testing " << mNumTestNodes << " nodes");
	for (int i = 0; i < node->mNumTestNodes; ++ i)
	{
		RenderableHierarchyNode *hnode = mTestNodes[node->mTestNodesIdx + i];
		const float dist = GetMinSquareDistance(hnode->GetBox());

		if (dist < minDist)
			minDist = dist;
	}

	distance = minDist;
#else
	distance = GetMinSquareDistance(node->GetBox());
#endif

	timeDistance.Exit();

	return distance;
}


float Bvh::GetMinSquareDistance(const BoundingBox &box) const
{
	float minDist = 1e25f;

	const Point3f *pts = (const Point3f *)box.getVertexData();
	const Point3f pos = mCamera->GetPosition();

	for (int i = 0; i < 8; ++ i)
	{
		const float newDist = DistanceSquared(pts[i], pos);

		if (newDist < minDist)
			minDist = newDist;
	}

	return minDist;
}


bool Bvh::ExportBinTree(const string &filename)
{
	ofstream stream(filename.c_str(), ifstream::binary);
	
	if (!stream.is_open()) return false;

	OUT1("exporting bvh");

	std::queue<BvhNode *> tStack;
	tStack.push(mRoot);

	while(!tStack.empty())
	{
		BvhNode *node = tStack.front();
		tStack.pop();
		
		if (node->IsLeaf())
		{
			//OUT1("l");
			ExportBinLeaf(stream, static_cast<BvhLeaf *>(node));
		}
		else
		{
			//OUT1("i");
			BvhInterior *interior = static_cast<BvhInterior *>(node);

			ExportBinInterior(stream, interior);
			
			tStack.push(interior->mFront);
			tStack.push(interior->mBack);
		}
	}

	OUT1("... finished");

	return true;
}


void Bvh::ExportBinLeaf(ofstream &stream, BvhLeaf *leaf)
{
	int type = TYPE_LEAF;
	int first = leaf->mFirst;
	int last = leaf->mLast;

	stream.write(reinterpret_cast<char *>(&type), sizeof(int));
	stream.write(reinterpret_cast<char *>(&first), sizeof(int));
	stream.write(reinterpret_cast<char *>(&last), sizeof(int));
}


void Bvh::ExportBinInterior(ofstream &stream, BvhInterior *interior)
{
	int type = TYPE_INTERIOR;
	stream.write(reinterpret_cast<char *>(&type), sizeof(int));
	stream.write(reinterpret_cast<char *>(&interior->mAxis), sizeof(char));
	stream.write(reinterpret_cast<char *>(&interior->mPosition), sizeof(float));
}


BvhLeaf *Bvh::ImportBinLeaf(ifstream &stream, BvhInterior *parent)
{
	BvhLeaf *leaf = new BvhLeaf(parent);

	int first, last;

	stream.read(reinterpret_cast<char *>(&first), sizeof(int));
	stream.read(reinterpret_cast<char *>(&last), sizeof(int));

	leaf->mFirst = first;
	leaf->mLast = last;

	return leaf;
}


BvhInterior *Bvh::ImportBinInterior(ifstream &stream, BvhInterior *parent)
{
	BvhInterior *interior = new BvhInterior(parent);

	stream.read(reinterpret_cast<char *>(&interior->mAxis), sizeof(char));
	stream.read(reinterpret_cast<char *>(&interior->mPosition), sizeof(float));

	return interior;
}


BvhNode *Bvh::LoadNextNode(ifstream &stream, BvhInterior *parent)
{
	int nodeType;
	stream.read(reinterpret_cast<char *>(&nodeType), sizeof(int));

	if (nodeType == TYPE_LEAF)
	{
		//OUT1("l");
		return ImportBinLeaf(stream, static_cast<BvhInterior *>(parent));
	}

	if (nodeType == TYPE_INTERIOR)
	{
		//OUT1("i");
		return ImportBinInterior(stream, static_cast<BvhInterior *>(parent));
	}

	return NULL;
}


Bvh *Bvh::LoadFromFile(const string &filename, const GeometryVector &geom)
{
	// export binary version of mesh
	queue<BvhNode *> tStack;
	ifstream stream(filename.c_str(), ifstream::binary);

	if (!stream.is_open()) return NULL;

	OUT1("loading bvh");
	Bvh *bvh = new Bvh();

	bvh->mGeometrySize = geom.size();
	bvh->mGeometry = new NodeGeometry*[bvh->mGeometrySize];

	for (size_t i = 0; i < bvh->mGeometrySize; ++ i)
		bvh->mGeometry[i] = geom[i];

	bvh->mRoot = bvh->LoadNextNode(stream, NULL);

	tStack.push(bvh->mRoot);
	bvh->mNumNodes = 1;

	while(!tStack.empty())
	{
		BvhNode *node = tStack.front();
		tStack.pop();

		if (!node->IsLeaf())
		{
			bvh->mNumNodes += 2;

			BvhInterior *interior = static_cast<BvhInterior *>(node);

            BvhNode *front = bvh->LoadNextNode(stream, interior);
			BvhNode *back = bvh->LoadNextNode(stream, interior);
	
			interior->mFront = front;
			interior->mBack = back;

			front->mDepth = interior->mDepth + 1;
			back->mDepth = interior->mDepth + 1;

			tStack.push(front);			
			tStack.push(back);
		}
	}

	OUT1("... finished loading " << bvh->mNumNodes << " nodes, updating boxes");

	//adjust bounding boxes
	bvh->UpdateBoxes(bvh->mRoot);
	bvh->UpdateNumLeaves(bvh->mRoot);

	// compute unique ids
	bvh->ComputeIds();
	// update the indices of the geometry so we can render interiors as well
	bvh->UpdateInteriors(bvh->mRoot);
	// do this once so at least the current node bounding box is tested
	bvh->RecomputeBounds();

	return bvh;
}


float Bvh::ComputeScreenSpaceProjection(BvhNode *node) const
{
#if 0
	int projection = 0;

	for (int i = 0; i < node->mNumTestNodes; ++ i)
	{
		RenderableHierarchyNode *n = mTestNodes[node->mTestNodesIdx + i];

		if(hnode->.mCache.mLastProjectionFrameId < mFrameId)
		{
			hnode->mCache.mProjection = 
				ComputeScreenSpaceProjection(hnode->GetBox());
			hnode->mCache.mLastProjectionFrameId = mCurrentFrameId;
		}
		
		projection += hnode->mCache.mProjection;
	}

	return projection;

#else

	if(node->mCache.mLastProjectionFrameId < mFrameId)
	{
		node->mCache.mProjection = 
			ComputeScreenSpaceProjection(node->GetBox());
		node->mCache.mLastProjectionFrameId = mFrameId;
	}

	return node->mCache.mProjection;
#endif
}


float Bvh::ComputeScreenSpaceProjection(const BoundingBox &box) const
{
#if 0
	const float dist = CalcDistance(box);
#else
	//const float dist = sqrt(GetMinSquareDistance(box));
	const float dist = Distance(mCamera->GetPosition(), box.getCenter());
#endif

	const float scale = dist ? dist * mScale : 1e-6;
	const float f = 1.0f / scale;
	const float coverage = f * f * box.getSurface() / 6.0f;

#if 0
	OUT1("scale: " << scale);
	OUT1("box: " << box.getSurface());
	OUT1("aspect: " << aspect);
	OUT1("fov: " << fov);
	OUT1("tan: " << tan(fov * mypi / 180.0f));
	//OUT1("w: " << w << " h: " << h);
	OUT1("cov: " << coverage << " pix: " << coverage * w * h << " f: " << f);
#endif
	return coverage;
}


void Bvh::RenderBoundingBox(BvhNode *node)
{
#if 1
	static BvhNodeContainer dummy(1);
	dummy[0] = node;
	RenderBoundingBoxes(dummy);
#else
	RenderBoxIndividual(node);
#endif
}


void Bvh::RenderBoundingBoxImmediate(const BoundingBox &box, const bool restartStrip)
{
	const Point3f *pU = box.getUpperPtr();
	const Point3f *pL = box.getLowerPtr();

	///////////
	//-- render AABB as triangle strips

	glVertex3f(pL->x, pL->y, pU->z);
	glVertex3f(pU->x, pL->y, pU->z);
	glVertex3f(pL->x, pU->y, pU->z);
	glVertex3f(pU->x, pU->y, pU->z);
	glVertex3f(pL->x, pU->y, pL->z); 
	glVertex3f(pU->x, pU->y, pL->z);
	glVertex3f(pL->x, pL->y, pL->z); 
	glVertex3f(pU->x, pL->y, pL->z);

	glPrimitiveRestartNV();

	//-- render second half of AABB
	glVertex3f(pL->x, pU->y, pU->z); 
	glVertex3f(pL->x, pU->y, pL->z);
	glVertex3f(pL->x, pL->y, pU->z);
	glVertex3f(pL->x, pL->y, pL->z);
	glVertex3f(pU->x, pL->y, pU->z); 
	glVertex3f(pU->x, pL->y, pL->z);
	glVertex3f(pU->x, pU->y, pU->z); 
	glVertex3f(pU->x, pU->y, pL->z);

	
	// restart for all but last node
	if (restartStrip)
	{
		glPrimitiveRestartNV();
	}
}


int Bvh::RenderBoundingBoxesImmediate(const BvhNodeContainer &nodes)
{
	timeSetup.Entry();
	int renderedBoxes = 0;
	
	glBegin(GL_TRIANGLE_STRIP);

	BvhNodeContainer::const_iterator bit, bit_end = nodes.end();

	for (bit = nodes.begin(); bit != bit_end; ++ bit)
	{
		BvhNode *node = *bit;

		renderedBoxes += node->mNumTestNodes;
#if 0
		const bool restartStrip = false;
		RenderBoundingBoxImmediate(node->GetBox(), restartStrip);
#else
		//OUT1("nodes: " << node->mTestNodes.size() << " l: " << node->IsLeaf());
		for (int i = 0; i < node->mNumTestNodes; ++ i)
		{
			RenderableHierarchyNode *testNode = mTestNodes[node->mTestNodesIdx + i];
			const BoundingBox &box = testNode->GetBox();

			// restart for all but last node
			const bool restartStrip = (bit != bit_end - 1) || (i != node->mNumTestNodes - 1);

			RenderBoundingBoxImmediate(box, restartStrip);
		}
#endif
	}
	glEnd();
	timeSetup.Exit();

	return renderedBoxes;
}


int Bvh::RenderBoundingBoxes(const BvhNodeContainer &nodes)
{
	// always render in immediate mode if there is only one box to render
	if (((nodes.size() == 1) && (nodes[0]->mNumTestNodes == 1)) || 
		!Settings::Global()->get_nvocc_use_index_arrays())
	{	
		// render bounding boxes individually
		return RenderBoundingBoxesImmediate(nodes);
	}
	else 
	{
		const int renderedBoxes = 
			PrepareBoundingBoxesWithDrawArrays(nodes);
		RenderBoundingBoxesWithDrawArrays(renderedBoxes);

		return renderedBoxes;
	}
}


int Bvh::PrepareBoundingBoxesWithDrawArrays(const BvhNodeContainer &nodes)
{
	//////
	//-- for the first time we come here ...

	if (!mIndices)
	{	// create list of indices
		CreateIndices();
	}

	if (sCurrentVboId == -1)
	{	
		// prepare the vbo 
		PrepareVertices();
	}

	///////////////

	timeSetup.Entry();

	int numNodes = 0;

	BvhNodeContainer::const_iterator nit, nit_end = nodes.end();

	for (nit = nodes.begin(); nit != nit_end; ++ nit)
	{
		BvhNode *node = *nit;

        const int numIndices = node->mNumTestNodes * sNumIndicesPerBox;
		
		// copy indices
		memcpy(mIndices + numNodes * sNumIndicesPerBox, 
			   mTestIndices + node->mIndicesPtr, 
#if 0 //ALIGN_INDICES
			   ((numIndices / 32) * 32 + 32) * sizeof(unsigned int));
#else
			   numIndices * sizeof(unsigned int));
#endif
		numNodes += node->mNumTestNodes;
	}

	timeSetup.Exit();

	return numNodes;
}


void Bvh::RenderBoundingBoxesWithDrawArrays(const int numNodes)
{
	//////
	//-- Rendering the vbo

	timeIssueDrawElements.Entry();

	if (useVbos)
		// set the vertex pointer to the vertex buffer
		glVertexPointer(3, GL_FLOAT, 0, (char *)NULL);		
	else
		glVertexPointer(3, GL_FLOAT, 0, mVertices);

	// we do use the last or the first index (they are generate and only used to connect strips)
	const int numElements = numNodes * sNumIndicesPerBox - 1;

	// don't render first degenerate index
	glDrawElements(GL_TRIANGLE_STRIP, numElements, GL_UNSIGNED_INT, mIndices + 1); 
	
	timeIssueDrawElements.Exit();
}


#define ALIGN_INDICES 1

void Bvh::CreateIndices()
{
	// collect bvh nodes
	BvhNodeContainer nodes;
	CollectNodes(mRoot, nodes);

	OUT1("creating new indices");

	int numMaxIndices = 0;

	BvhNodeContainer::const_iterator lit, lit_end = nodes.end();

	for (lit = nodes.begin(); lit != lit_end; ++ lit)
	{
		int offset = (*lit)->mNumTestNodes * sNumIndicesPerBox;
#if ALIGN_INDICES
		// align with 32
		offset = (offset / 32) * 32 + 32; 
#endif
		numMaxIndices += offset;
	}


	OUT1("creating global indices buffer");

	if (mIndices) delete [] mIndices;
	if (mTestIndices) delete [] mTestIndices;

	// global buffer: create it once so we don't have
	// to allocate memory from the chunks of the node
	mIndices = new unsigned int[numMaxIndices];

	// create new index buffer for the individual nodes
	mTestIndices = new unsigned int[numMaxIndices];
	
	mCurrentIndicesPtr = 0;

	for (lit = nodes.begin(); lit != lit_end; ++ lit)
	{
		BvhNode *node = *lit;
		
		// resize array
		node->mIndicesPtr = mCurrentIndicesPtr;

		int numIndices = 0;

		// the bounding box of the test nodes are rendered instead of the root node
		// => store their indices
		for (int i = 0; i < node->mNumTestNodes; ++ i, numIndices += sNumIndicesPerBox)
		{
			RenderableHierarchyNode *testNode = mTestNodes[node->mTestNodesIdx + i];

			// add indices to root node
			for (int j = 0; j < sNumIndicesPerBox; ++ j)
			{
				mTestIndices[mCurrentIndicesPtr + numIndices + j] = sIndices[j] + testNode->GetId() * 8;
			}
		}

		// align with 32
#if ALIGN_INDICES
		const int offset = (numIndices / 32) * 32 + 32;
#else
		const int offset = numIndices;
#endif
		mCurrentIndicesPtr += offset;
	}
}


void Bvh::ComputeIds()
{
	// collect all nodes, also the nodes from local bvh
	// warning: root nodes local bvh are not in there, as they 
	// are equivalent geometry bvh leaves
	HierarchyNodeContainer nodes;
	CollectAllNodes(mRoot, nodes);

	// assign ids to all nodes of the "regular" hierarchy
	int i = 0;
	HierarchyNodeContainer::const_iterator lit, lit_end = nodes.end();

	for (lit = nodes.begin(); lit != lit_end; ++ lit, ++ i)
	{
		(*lit)->SetId(i);
	}
}


void Bvh::PrepareVertices()
{
	// collect all nodes, also the nodes from local bvh
	// warning: root nodes local bvh are not in there, as they 
	// are equivalent geometry bvh leaves
	HierarchyNodeContainer nodes;

	nodes.reserve(GetNumNodes());
	CollectAllNodes(mRoot, nodes);

	const unsigned int bufferSize = 8 * (int)nodes.size();
	mVertices = new Point3f[bufferSize];
	
	int i = 0;
	
	// store bounding box vertices
	HierarchyNodeContainer::const_iterator lit, lit_end = nodes.end();

	for (lit = nodes.begin(); lit != lit_end; ++ lit, i += 8)
	{
		RenderableHierarchyNode *node = *lit;
		const Point3f *vertices = (const Point3f *)node->GetBox().getVertexData();

		for (int j = 0; j < 8; ++ j)
		{	
			((Point3f *)mVertices)[node->GetId() * 8 + j] = vertices[j];
		}
	}

	if (useVbos)
	{
		glGenBuffersARB(1, &sCurrentVboId);
		glBindBufferARB(GL_ARRAY_BUFFER_ARB, sCurrentVboId);
		glVertexPointer(3, GL_FLOAT, 0, (char *)NULL);
	
		glBufferDataARB(GL_ARRAY_BUFFER_ARB, 
			            bufferSize * sizeof(Point3f), 
			            mVertices, 
						GL_STATIC_DRAW_ARB);

		//glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);

		// data handled by graphics driver from now on
		DELAPTR(mVertices);

		OUT1("***** created vbos *********");
	}
	else
	{
		sCurrentVboId = 0;
        glVertexPointer(3, GL_FLOAT, 0, mVertices);

		OUT1("created vertices");
	}
}


void Bvh::SetMaxDepthForTestingChildren(const int maxDepth) 
{
	if (maxDepth != mMaxDepthForTestingChildren)
	{
		mMaxDepthForTestingChildren = maxDepth;
		RecomputeBounds();
	}
}


void Bvh::SetAreaRatioThresholdForTestingChildren(const float ratio) 
{
	if (ratio != mAreaRatioThreshold)
	{
		mAreaRatioThreshold = ratio;
		RecomputeBounds();
	}
}


void Bvh::SetVolRatioThresholdForTestingChildren(const float ratio) 
{
	if (ratio != mVolRatioThreshold)
	{
		mVolRatioThreshold = ratio;
		RecomputeBounds();
	}
}


void Bvh::SetUseTighterBoundsForTests(const bool tighterBoundsForTests)
{
	if (mUseTighterBoundsOnlyForLeafTests != tighterBoundsForTests)
	{
		mUseTighterBoundsOnlyForLeafTests = tighterBoundsForTests;
		RecomputeBounds();
	}
}


void Bvh::SetCollectTighterBoundsWithMaxLevel(bool t)
{
	if (mCollectTighterBoundsWithMaxLevel != t)
	{
		mCollectTighterBoundsWithMaxLevel = t;
		RecomputeBounds();
	}
}


void Bvh::RecomputeBounds()
{
	// clear old list
	mTestNodes.clear();

	// collect all nodes
	BvhNodeContainer nodes;
	CollectNodes(mRoot, nodes);

	OUT1("recomputing bounds, children will be tested in depth " << mMaxDepthForTestingChildren);

	if (mCollectTighterBoundsWithMaxLevel)
		OUT1("creating tighter bounds using max level");
	else
		OUT1("creating tighter bounds using best set");

	int success = 0;
	BvhNodeContainer::const_iterator lit, lit_end = nodes.end();

	for (lit = nodes.begin(); lit != lit_end; ++ lit)
	{
		BvhNode *node = *lit;

		if (mCollectTighterBoundsWithMaxLevel)
		{
			//OUT1("creating tighter bounds using max level");

			// recreate list of nodes that will be tested as a proxy ...
			if (CreateNodeRenderList(node))
				++ success;
		}
		else
		{
			//OUT1("creating tighter bounds using best set");

			HierarchyNodeContainer hnodes;
			CollectBestNodes(node, hnodes);	

			// the new test nodes are added at the end of the vector
			node->mTestNodesIdx = mTestNodes.size();

			HierarchyNodeContainer::const_iterator cit;

			// use the found nodes as nodes during the occlusion tests
			for (cit = hnodes.begin(); cit != hnodes.end(); ++ cit)
			{
				RenderableHierarchyNode *child = *cit;
				mTestNodes.push_back(child);
			}

			node->mNumTestNodes = (int)hnodes.size();
		}
		//OUT1("testnodes: " << node->mNumTestNodes);
	}

	const float p = 100.0f * (float)success / nodes.size();

	OUT1("created tighter bounds for " << p << " percent of the nodes");

	// recreate indices used for indirect mode rendering
	if (mIndices)
	{
		CreateIndices();
	}
}

	
bool Bvh::CreateNodeRenderList(BvhNode *node)
{
	HierarchyNodeContainer children;

	if (mUseTighterBoundsOnlyForLeafTests && !node->IsLeaf())
	{
		children.push_back(node);
	}
	else
	{
		// collect nodes that will be tested instead of the leaf node
		// in order to get a tighter bounding box test
		CollectNodes(node, mMaxDepthForTestingChildren, children);
	}


	// using the tighter bounds is not feasable in case
	// that the tighter bounds represent nearly the same projected area
	// as the old bounding box. Find this out using either volume or area
	// heuristics

	float vol = 0;
	float area = 0;

	HierarchyNodeContainer::const_iterator cit;

	for (cit = children.begin(); cit != children.end(); ++ cit)
	{
		RenderableHierarchyNode *c = *cit;

		area += c->GetBox().getSurface();
		vol += c->GetBox().getVolume();
	}

	const float volRatio = vol / node->GetBox().getVolume();
	const float areaRatio = area / node->GetBox().getSurface();
	
	bool success;

	if ((areaRatio < mAreaRatioThreshold) &&
		(volRatio < mVolRatioThreshold))
	{
		success = true;
	}
	else
	{
		// hack: only store node itself
		children.clear();
		children.push_back(node);

		success = false;
	}

	// the new test nodes are added at the end of the vector
	node->mTestNodesIdx = mTestNodes.size();

	// use the found nodes as nodes during the occlusion tests
	for (cit = children.begin(); cit != children.end(); ++ cit)
	{
		RenderableHierarchyNode *child = *cit;
		mTestNodes.push_back(child);
	}

	node->mNumTestNodes = (int)children.size();

	return success;
}


void Bvh::ResetNodeClassifications()
{
	BvhNodeContainer nodes;

	nodes.reserve(GetNumNodes());
	CollectNodes(mRoot, nodes);

	BvhNodeContainer::const_iterator lit, lit_end = nodes.end();

	for (lit = nodes.begin(); lit != lit_end; ++ lit)
	{
		(*lit)->ResetVisibility();
	}
}


int Bvh::PostProcessLeaves(BvhLeafContainer &leaves)
{
	OUT1("creating tighter bounds for leaves");

	BvhLeafContainer::const_iterator lit, lit_end = leaves.end();

	int i = 0;
	int tighter = 0;

	for (lit = leaves.begin(); lit != lit_end; ++ lit, ++ i)
	{
		BvhLeaf *leaf = *lit;

		int triangleCount;
		// collect the triangles from the stored geometry
		Point3f *triangles = CollectTriangles(leaf, triangleCount);

		if (CreateTighterBoundsForLeaf(leaf, triangles, triangleCount))
			++ tighter;

		leaf->mArea = ComputeGeometryArea(leaf, triangles, triangleCount);

		delete [] triangles;

		if (i % 1000 == 999)
			OUT1(i << " leaves processed ");
	}

	const float p = 100.0f * tighter / leaves.size();
	OUT1("created tighter bounds for " << p << " percent of the leaves");

	return tighter;
}


bool Bvh::CreateTighterBoundsForLeaf(BvhLeaf *leaf, Point3f *triangles, const int triangleCount)
{
	// create a local bvh over the triangles
	if (leaf->mTriangleBvh) delete leaf->mTriangleBvh;

	leaf->mTriangleBvh = new TriangleBvh((Triangle *)triangles, triangleCount);

	const int maxDepth = Settings::Global()->get_nvocc_local_bvh_max_depth();
	const int maxTriangles = Settings::Global()->get_nvocc_local_bvh_max_triangles();
	const int minArea = Settings::Global()->get_nvocc_local_bvh_min_area();
	const int splitType = Settings::Global()->get_nvocc_local_bvh_split_type();

	leaf->mTriangleBvh->SetMaxDepth(maxDepth);
	leaf->mTriangleBvh->SetMaxTriangles(maxTriangles);
	leaf->mTriangleBvh->SetMinAreaRatio(minArea);
	leaf->mTriangleBvh->SetSplitType(splitType);
	
	leaf->mTriangleBvh->SetVerboseOutput(false);

	// construct
	leaf->mTriangleBvh->Construct();

	return true;
}


float Bvh::ComputeGeometryArea(BvhLeaf *leaf, Point3f *triangles, const int triangleCount) const
{
	float area = 0;

	for (int i = 0; i < triangleCount; ++ i)
	{
		area += ((Triangle *)triangles)[i].GetArea();
	}

	return area;
}


Point3f *Bvh::CollectTriangles(BvhLeaf *leaf, int &triangleCount)
{
	// count triangles in the leaf
	triangleCount = 0;

	for (int geomIdx = leaf->mFirst; geomIdx <= leaf->mLast; ++ geomIdx)
	{
		Shape3D *sh = static_cast<Shape3D *>(mGeometry[geomIdx]->GetNode());
		IndexedTriangleStripArray *strips = static_cast<IndexedTriangleStripArray *>(sh->getGeometry());

		triangleCount += strips->getNumTriangles();
	}

	Triangle *triangles = new Triangle[triangleCount];
	//OUT1("The leaf contains " << triangleCount << " triangles in " << (int)geometry.size() << " nodes");

	int currentIdx = 0;
	
	for (int geomIdx = leaf->mFirst; geomIdx <= leaf->mLast; ++ geomIdx)
	{
		NodeGeometry *geom = mGeometry[geomIdx];

		Shape3D *sh = static_cast<Shape3D *>(geom->GetNode());
		
		IndexedTriangleStripArray *strips = 
			static_cast<IndexedTriangleStripArray *>(sh->getGeometry());

		// get indices of containted triangles (strips are converted before)
		int *indices;
	
		int tCount = strips->getNumTriangles();
		int indexCount = tCount * 3;
		
		// compute triangles from strips
		Point3f *coordinates = (Point3f *)strips->getCoordRef3f();
		strips->getTriangleCoordinateIndices(0, &indices, indexCount);

		for (int i = 0; i < tCount; ++ i)
		{
			triangles[i + currentIdx].mVertices[0] = coordinates[indices[i * 3 + 0]]; 
			triangles[i + currentIdx].mVertices[1] = coordinates[indices[i * 3 + 1]]; 
			triangles[i + currentIdx].mVertices[2] = coordinates[indices[i * 3 + 2]]; 

			if (geom->mMatId != NV_RenderPredictorRenderAction::NO_MAT)
			{
				const Matrix4f &m = mRenderer->app_mats[geom->mMatId];
				Transform3D tf(m);
				tf.transform(triangles[i + currentIdx].mVertices[0]);
				tf.transform(triangles[i + currentIdx].mVertices[1]);
				tf.transform(triangles[i + currentIdx].mVertices[2]);
			}
		}

		currentIdx += tCount;
	}

	return (Point3f *)triangles;
}


int Bvh::CountTriangles(BvhLeaf *leaf) const
{
	int triangleCount = 0;
	
	for (int i = leaf->mFirst; i <= leaf->mLast; ++ i)
	{
		triangleCount += mGeometry[i]->mNumTriangles;
	}

	return triangleCount;
}


void Bvh::ComputeBvhStats() 
{
	std::stack<BvhNode *> nodeStack;
	nodeStack.push(mRoot);

	while (!nodeStack.empty()) 
	{
		BvhNode *node = nodeStack.top();
		nodeStack.pop();

		if (node->IsLeaf()) 
		{
			BvhLeaf *leaf = static_cast<BvhLeaf *>(node);

			mBvhStats.mLeafSA += leaf->mBox.getSurface();
			mBvhStats.mLeafVol += leaf->mBox.getVolume();

			TriangleBvh::TriangleBvhStats tStats;
			leaf->mTriangleBvh->GetBvhStats(tStats);

			mBvhStats.mBoundsLeafSA += tStats.mLeafSA;
			mBvhStats.mBoundsInteriorSA += tStats.mInteriorSA;

			mBvhStats.mBoundsLeafVol += tStats.mLeafVol;
			mBvhStats.mBoundsInteriorVol += tStats.mInteriorVol;

			mBvhStats.mBoundsLeavesCount += leaf->mTriangleBvh->GetNumLeaves();
		} 
		else 
		{
			mBvhStats.mInteriorSA += node->mBox.getSurface();
			mBvhStats.mInteriorVol += node->mBox.getVolume();

			BvhInterior *interior = static_cast<BvhInterior *>(node);
			
			nodeStack.push(interior->mBack);
			nodeStack.push(interior->mFront);
		}
	}

	mBvhStats.mGeometryRatio = mGeometrySize / (float)GetNumLeaves();
	mBvhStats.mTriangleRatio = mBvhStats.mTriangles / (float)GetNumLeaves();
}


void Bvh::PrintBvhStats() const
{
	//OUT1("triangle bvh surface " << mBvhStats.mBoundsLeafSA);
	//OUT1("triangle bvh volume " << mBvhStats.mBoundsLeafVol);

	OUT1("bvh stats:");
	OUT1("interiorNodesSA = " << mBvhStats.mInteriorSA / mRoot->mBox.getSurface());
	OUT1("leafNodesSA = " << mBvhStats.mLeafSA / mRoot->mBox.getSurface());
	OUT1("interiorNodesVolume = " << mBvhStats.mInteriorVol / mRoot->mBox.getVolume());
	OUT1("leafNodesVolume = " << mBvhStats.mLeafVol / mRoot->mBox.getVolume() << "\n");

	OUT1("boundsInteriorNodesSA = " << mBvhStats.mBoundsInteriorSA / mBvhStats.mLeafSA);
	OUT1("boundsLeafNodesSA = " << mBvhStats.mBoundsLeafSA / mBvhStats.mLeafSA);
	OUT1("boundsInteriorNodesVolume = " << mBvhStats.mBoundsInteriorVol / mBvhStats.mLeafVol);
	OUT1("boundsLeafNodesVolume = " << mBvhStats.mBoundsLeafVol / mBvhStats.mLeafVol << "\n");

	OUT1("boundsLeaves: " <<  (float)mBvhStats.mBoundsLeavesCount / (float)GetNumLeaves() << "\n");
	OUT1("geometry per leaf: " <<  mBvhStats.mGeometryRatio);
	OUT1("triangles per leaf: " <<  mBvhStats.mTriangleRatio);
	OUT1("");
}


static void RenderBoxForViz(const BoundingBox &box)
{
	glBegin(GL_LINE_LOOP);
	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
	glEnd();

	glBegin(GL_LINE_LOOP);
	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
	glEnd();

	glBegin(GL_LINE_LOOP);
	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
	glEnd();

	glBegin(GL_LINE_LOOP);
	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
	glEnd();

	glBegin(GL_LINE_LOOP);
	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
	glEnd();

	glBegin(GL_LINE_LOOP);
	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);

	glEnd();
}


void Bvh::RenderBoundingBoxesForViz(const int mode)
{
	BvhLeafContainer leaves;
	leaves.reserve(mRoot->GetNumLeaves());

	CollectLeaves(mRoot, leaves);

	BvhLeafContainer::const_iterator it, it_end = leaves.end();

	for (it = leaves.begin(); it != it_end; ++ it)
	{
		BvhLeaf *leaf = *it;

		if ((mode == 0) || (mode == 2))
		{
			glColor3f(1.0f, 1.0f, 1.0f);
            RenderBoxForViz(leaf->GetBox());
		}
		
		if ((mode == 1) || (mode == 2)) 
		{
			glColor3f(1.0f, 0, 0);
			for (int i = 0; i < leaf->mNumTestNodes; ++ i)
			{
				RenderBoxForViz(mTestNodes[leaf->mTestNodesIdx + i]->GetBox());
			}
		}
	}
}


int Bvh::CountTriangles(BvhNode *node) const
{
	int numTriangles = 0;

	for (int i = node->mFirst; i <= node->mLast; ++ i)
	{
		numTriangles += mGeometry[i]->mNumTriangles;
	}

	return numTriangles;
}


float Bvh::GetGeometryArea(BvhNode *node) const
{
	return node->mArea;
}


float Bvh::GetAvgDepth() const
{
	return mAvgDepth;
}


static float GetNodePriority(RenderableHierarchyNode *node)
{
	if (node->IsLeaf())
		return 0.0f;

	float result;

	if (node->GetType() == BVH_NODE)
	{
		BvhInterior *interior = static_cast<BvhInterior *>(node);
		
		// evaluate the priority of this node
		const float correctionFactor = 0.0f;
		//  1.0f*interior->box.FaceArea(node->axis);

		result = 
			interior->GetBox().getSurface() - 
			(interior->GetBack()->GetBox().getSurface() + interior->GetFront()->GetBox().getSurface()) + 
			correctionFactor;
	}
	else
	{
		TriangleBvhInterior *interior = static_cast<TriangleBvhInterior *>(node);
		
		// evaluate the priority of this node
		const float correctionFactor = 0.0f;
		//  1.0f*interior->box.FaceArea(node->axis);

		result = 
			interior->GetBox().getSurface() - 
			(interior->GetBack()->GetBox().getSurface() + interior->GetFront()->GetBox().getSurface()) + 
			correctionFactor;
	}

	return result;
}


struct NodeEntry
{
	NodeEntry(RenderableHierarchyNode *node, float p): mNode(node), mPriority(p)
	{}

	RenderableHierarchyNode *mNode;
	float mPriority;
};


bool operator<(const NodeEntry &a, const NodeEntry &b)
{
	return a.mPriority < b.mPriority;
}


void Bvh::CollectBestNodes(RenderableHierarchyNode *node, HierarchyNodeContainer &nodes)
{
	int maxNodes = 32;
	priority_queue<NodeEntry> nodeStack;

	nodeStack.push(NodeEntry(node, GetNodePriority(node)));
	nodes.clear();

	float SA = node->GetBox().getSurface();
	float saThreshold = 1.1f * SA;
	int numNodes = 1;
	OUT1(numNodes << " " << SA);

	while (!nodeStack.empty() &&
		int(nodes.size() + nodeStack.size()) < maxNodes)
	{
		NodeEntry entry = nodeStack.top();
		nodeStack.pop();

		RenderableHierarchyNode *current = entry.mNode;

		if (current->IsLeaf())
		{
			// check if there further subdivision on triangle level
			if ((current->GetType() == BVH_NODE) && (((BvhLeaf *)current)->mTriangleBvh->GetNumNodes() > 1))
			{
				// push back the root of the local bvh
				nodeStack.push(NodeEntry(((BvhLeaf *)current)->mTriangleBvh->GetRoot(), entry.mPriority));
			}
			else // terminate traversal
			{
				nodes.push_back(current);
			}
		}
		else // interior node
		{
			// surface area of child nodes
			SA -= entry.mPriority;
		
			// finish the search if we have too much fillrate increase
			if (SA > saThreshold)
			{
				OUT1("break at " << SA << " > " << saThreshold);
				nodes.push_back(current);
			}
			else
			{
				++ numNodes;

				OUT1(numNodes << " " << SA << " " << entry.mPriority);

				if (current->GetType() == BVH_NODE)
				{
					BvhInterior *interior = static_cast<BvhInterior *>(current);

					nodeStack.push(NodeEntry(interior->GetBack(), 
						GetNodePriority(interior->GetBack())));
					nodeStack.push(NodeEntry(interior->GetFront(),
						GetNodePriority(interior->GetFront())));
				}
				else
				{
					TriangleBvhInterior *interior = static_cast<TriangleBvhInterior *>(current);

					nodeStack.push(NodeEntry(interior->GetBack(), 
						GetNodePriority(interior->GetBack())));
					nodeStack.push(NodeEntry(interior->GetFront(),
						GetNodePriority(interior->GetFront())));
				}
			}
		}
	}

	while (!nodeStack.empty()) 
	{
		NodeEntry entry = nodeStack.top();
	
		nodeStack.pop();
		RenderableHierarchyNode *current = entry.mNode;
		nodes.push_back(current);
	}
} 

#endif