#include <stack>
#include <time.h>
#include <iomanip>

#include "BvHierarchy.h"
#include "ViewCell.h"
#include "Plane3.h"
#include "Mesh.h"
#include "common.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include "Exporter.h"
#include "Plane3.h"
#include "ViewCellsManager.h"
#include "Beam.h"
#include "VspTree.h"
#include "HierarchyManager.h"


namespace GtpVisibilityPreprocessor {


#define PROBABILIY_IS_BV_VOLUME 1
#define USE_FIXEDPOINT_T 0
#define USE_VOLUMES_FOR_HEURISTICS 1
#define TEST_POWERPLANT 0
#define USE_BETTER_RENDERCOST_EST 0

//int BvhNode::sMailId = 10000;
//int BvhNode::sReservedMailboxes = 1;

BvHierarchy *BvHierarchy::BvhSubdivisionCandidate::sBvHierarchy = NULL;


/// sorting operator
inline static bool ilt(Intersectable *obj1, Intersectable *obj2)
{
	return obj1->mId < obj2->mId;
}


/// sorting operator
inline static bool smallerSize(Intersectable *obj1, Intersectable *obj2)
{
	return obj1->GetBox().SurfaceArea() < obj2->GetBox().SurfaceArea();
}



/***************************************************************/
/*              class BvhNode implementation                   */
/***************************************************************/

BvhNode::BvhNode(): 
mParent(NULL), 
mTimeStamp(0),
mRenderCost(-1)

{
	
}

BvhNode::BvhNode(const AxisAlignedBox3 &bbox):
mParent(NULL), 
mBoundingBox(bbox), 
mTimeStamp(0),
mRenderCost(-1)
{
}


BvhNode::BvhNode(const AxisAlignedBox3 &bbox, BvhInterior *parent):
mBoundingBox(bbox), 
mParent(parent), 
mTimeStamp(0),
mRenderCost(-1)
{
}


bool BvhNode::IsRoot() const
{
	return mParent == NULL;
}


BvhInterior *BvhNode::GetParent()
{
	return mParent;
}


void BvhNode::SetParent(BvhInterior *parent)
{
	mParent = parent;
}


int BvhNode::GetRandomEdgePoint(Vector3 &point,
								Vector3 &normal)
{
	// get random edge
	const int idx = Random(12); 
	Vector3 a, b;
	mBoundingBox.GetEdge(idx, &a, &b);
	
	const float w = RandomValue(0.0f, 1.0f);

	point = a * w + b * (1.0f - w);

	// TODO
	normal = Vector3(0);

	return idx;
}



/******************************************************************/
/*              class BvhInterior implementation                  */
/******************************************************************/


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox):
BvhNode(bbox), 
mSubdivisionCandidate(NULL),
mGlList(0)
{
  mActiveNode = this;
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent):
  BvhNode(bbox, parent),
  mGlList(0)
  
{
	mActiveNode = this;
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, 
				 BvhInterior *parent, 
				 const int numObjects):
  BvhNode(bbox, parent),
  mGlList(0)

{
	mObjects.reserve(numObjects);
	mActiveNode = this;
}


bool BvhLeaf::IsLeaf() const 
{ 
	return true; 
}


BvhLeaf::~BvhLeaf() 
{ 
}


void BvhLeaf::CollectObjects(ObjectContainer &objects)
{
	ObjectContainer::const_iterator oit, oit_end = mObjects.end();
	for (oit = mObjects.begin(); oit != oit_end; ++ oit)
	{
		objects.push_back(*oit);
	}
}



/******************************************************************/
/*              class BvhInterior implementation                  */
/******************************************************************/


BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox):
BvhNode(bbox), mFront(NULL), mBack(NULL)
{
}


BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox, BvhInterior *parent):
BvhNode(bbox, parent), mFront(NULL), mBack(NULL)
{
}


void BvhInterior::ReplaceChildLink(BvhNode *oldChild, BvhNode *newChild) 
{
	if (mBack == oldChild)
		mBack = newChild;
	else
		mFront = newChild;
}


bool BvhInterior::IsLeaf() const 
{ 
	return false; 
}


BvhInterior::~BvhInterior() 
{
	DEL_PTR(mFront); 
	DEL_PTR(mBack);
}


void BvhInterior::SetupChildLinks(BvhNode *front, BvhNode *back) 
{
    mBack = back;
    mFront = front;
}


void BvhInterior::CollectObjects(ObjectContainer &objects)
{
	mFront->CollectObjects(objects);
	mBack->CollectObjects(objects);
}


/*******************************************************************/
/*                  class BvHierarchy implementation               */
/*******************************************************************/


BvHierarchy::BvHierarchy():
mRoot(NULL),
mTimeStamp(1),
mIsInitialSubdivision(false)
{
	ReadEnvironment();
	mSubdivisionCandidates = new SortableEntryContainer;
//	for (int i = 0; i < 4; ++ i)
//		mSortedObjects[i] = NULL;
}


BvHierarchy::~BvHierarchy()
{
	// delete the local subdivision candidates
	DEL_PTR(mSubdivisionCandidates);

	// delete the presorted objects
	/*for (int i = 0; i < 4; ++ i)
	{
		DEL_PTR(mSortedObjects[i]);
	}*/
	
	// delete the tree
	DEL_PTR(mRoot);
}


void BvHierarchy::ReadEnvironment()
{
	bool randomize = false;
	Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.randomize", randomize);
	 
	// initialise random generator for heuristics
	if (randomize)
		Randomize();

	//////////////////////////////
	//-- termination criteria for autopartition

	Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.maxLeaves", mTermMaxLeaves);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.minObjects", mTermMinObjects);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.minRays", mTermMinRays);
	Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.minProbability", mTermMinProbability);
    Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.missTolerance", mTermMissTolerance);


	//////////////////////////////
	//-- max cost ratio for early tree termination

	Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.maxCostRatio", mTermMaxCostRatio);
	Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.minGlobalCostRatio",
		mTermMinGlobalCostRatio);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.globalCostMissTolerance", 
		mTermGlobalCostMissTolerance);


	//////////////////////////////
	//-- factors for subdivision heuristics

	// if only the driving axis is used for splits
	Environment::GetSingleton()->GetBoolValue("BvHierarchy.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
	Environment::GetSingleton()->GetFloatValue("BvHierarchy.maxStaticMemory", mMaxMemory);
	Environment::GetSingleton()->GetBoolValue("BvHierarchy.useCostHeuristics", mUseCostHeuristics);
	Environment::GetSingleton()->GetBoolValue("BvHierarchy.useSah", mUseSah);

    char subdivisionStatsLog[100];
	Environment::GetSingleton()->GetStringValue("BvHierarchy.subdivisionStats", subdivisionStatsLog);
	mSubdivisionStats.open(subdivisionStatsLog);

	Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
	Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.useGlobalSorting", mUseGlobalSorting);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.minRaysForVisibility", mMinRaysForVisibility);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.maxTests", mMaxTests);
	Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.useInitialSubdivision", mApplyInitialPartition);
	Environment::GetSingleton()->GetIntValue("BvHierarchy.Construction.Initial.minObjects", mInitialMinObjects);
	Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.Initial.maxAreaRatio", mInitialMaxAreaRatio);
	Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.Initial.minArea", mInitialMinArea);

	//mMemoryConst = (float)(sizeof(VspLeaf) + sizeof(VspViewCell));
	//mMemoryConst = (float)sizeof(BvhLeaf);
	mMemoryConst = 16;//(float)sizeof(ObjectContainer);

    mUseBboxAreaForSah = true;

	/////////////
	//-- debug output

	Debug << "******* Bvh hierarchy options ******** " << endl;
    Debug << "max depth: " << mTermMaxDepth << endl;
	Debug << "min probabiliy: " << mTermMinProbability<< endl;
	Debug << "min objects: " << mTermMinObjects << endl;
	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	Debug << "miss tolerance: " << mTermMissTolerance << endl;
	Debug << "max leaves: " << mTermMaxLeaves << endl;
	Debug << "randomize: " << randomize << endl;
	Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
	Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
	Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
	Debug << "max memory: " << mMaxMemory << endl;
	Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
	Debug << "use surface area heuristics: " << mUseSah << endl;
	Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
	//Debug << "split borders: " << mSplitBorder << endl;
	Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
	Debug << "use global sort: " << mUseGlobalSorting << endl;
	Debug << "minimal rays for visibility: " << mMinRaysForVisibility << endl;
	Debug << "bvh mem const: " << mMemoryConst << endl;
	Debug << "apply initial partition: " << mApplyInitialPartition << endl;
	Debug << "min objects: " << mInitialMinObjects << endl;
	Debug << "max area ratio: " << mInitialMaxAreaRatio << endl;
	Debug << "min area: " << mInitialMinArea << endl;

	Debug << endl;
}


void BvHierarchy::AssociateObjectsWithLeaf(BvhLeaf *leaf)
{
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		(*oit)->mBvhLeaf = leaf;
	}
}


static int CountRays(const ObjectContainer &objects)
{
	int nRays = 0;

	ObjectContainer::const_iterator oit, oit_end = objects.end();

	// warning: not exact number (there can be rays counted twice)
	// otherwise we would have to traverse trough all rays
	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		nRays += (int)(*oit)->GetOrCreateRays()->size();
	}

	return nRays;
}

									
float BvHierarchy::GetViewSpaceVolume() const
{
	return mViewCellsManager->GetViewSpaceBox().GetVolume();
}


BvhInterior *BvHierarchy::SubdivideNode(const BvhSubdivisionCandidate &sc,
										BvhTraversalData &frontData,
										BvhTraversalData &backData)
{ 
	// fill view cells cache
	mNodeTimer.Entry();

	const BvhTraversalData &tData = sc.mParentData;
	BvhLeaf *leaf = tData.mNode;

	AxisAlignedBox3 parentBox = leaf->GetBoundingBox();

	// update stats: we have two new leaves
	mBvhStats.nodes += 2;

	if (tData.mDepth > mBvhStats.maxDepth)
	{
		mBvhStats.maxDepth = tData.mDepth;
	}

	// add the new nodes to the tree
	BvhInterior *node = new BvhInterior(parentBox, leaf->GetParent());


	//////////////////
	//-- create front and back leaf

	AxisAlignedBox3 fbox = EvalBoundingBox(sc.mFrontObjects, &parentBox);
	AxisAlignedBox3 bbox = EvalBoundingBox(sc.mBackObjects, &parentBox);

	BvhLeaf *back = new BvhLeaf(bbox, node, (int)sc.mBackObjects.size());
	BvhLeaf *front = new BvhLeaf(fbox, node, (int)sc.mFrontObjects.size());

	BvhInterior *parent = leaf->GetParent();

	// replace a link from node's parent
	if (parent)
	{
		parent->ReplaceChildLink(leaf, node);
		node->SetParent(parent);
	}
	else // no parent => this node is the root
	{
		mRoot = node;
	}

	// and setup child links
	node->SetupChildLinks(front, back);

	++ mBvhStats.splits;

 
	////////////////////////////////////////
	//-- fill front and back traversal data with the new values

	frontData.mDepth = backData.mDepth = tData.mDepth + 1;

	frontData.mNode = front;
	backData.mNode = back;

	backData.mSampledObjects = new ObjectContainer();
	frontData.mSampledObjects = new ObjectContainer();

	*backData.mSampledObjects = sc.mSampledBackObjects;
	*frontData.mSampledObjects = sc.mSampledFrontObjects;

	back->mObjects = sc.mBackObjects;
	front->mObjects = sc.mFrontObjects;

	// if the number of rays is too low, no assumptions can be made
	// (=> switch to surface area heuristics?)
	frontData.mNumRays = CountRays(sc.mSampledFrontObjects);
	backData.mNumRays = CountRays(sc.mSampledBackObjects);

	AssociateObjectsWithLeaf(back);
	AssociateObjectsWithLeaf(front);
  
	////////////
	//-- compute pvs correction to cope with undersampling

	frontData.mPvs = (float)sc.mNumFrontViewCells;//(float)CountViewCells(sc.mSampledFrontObjects);
	backData.mPvs = (float)sc.mNumBackViewCells;//(float)CountViewCells(sc.mSampledBackObjects);

	frontData.mCorrectedPvs = sc.mCorrectedFrontPvs;
	backData.mCorrectedPvs = sc.mCorrectedBackPvs;


	// compute probability of this node being visible, 
	// i.e., volume of the view cells that can see this node
	frontData.mVolume = sc.mVolumeFrontViewCells;//EvalViewCellsVolume(sc.mSampledFrontObjects) / GetViewSpaceVolume();
	backData.mVolume = sc.mVolumeBackViewCells;//EvalViewCellsVolume(sc.mSampledBackObjects) / GetViewSpaceVolume();

	frontData.mCorrectedVolume = sc.mCorrectedFrontVolume;
	backData.mCorrectedVolume = sc.mCorrectedBackVolume;


    // how often was max cost ratio missed in this branch?
	frontData.mMaxCostMisses = sc.GetMaxCostMisses();
	backData.mMaxCostMisses = sc.GetMaxCostMisses();

	// set the time stamp so the order of traversal can be reconstructed
	node->SetTimeStamp(mHierarchyManager->mTimeStamp ++);
	 
	// assign the objects in sorted order
	if (mUseGlobalSorting)
	{
		AssignSortedObjects(sc, frontData, backData);
	}

	mNodeTimer.Exit();

	// return the new interior node
	return node;
}


BvhNode *BvHierarchy::Subdivide(SplitQueue &tQueue,
								SubdivisionCandidate *splitCandidate,
								const bool globalCriteriaMet
								,vector<SubdivisionCandidate *> &dirtyList
								)
{
	mSubdivTimer.Entry();

	BvhSubdivisionCandidate *sc = 
		static_cast<BvhSubdivisionCandidate *>(splitCandidate);
	BvhTraversalData &tData = sc->mParentData;

	BvhNode *currentNode = tData.mNode;

	if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
	{	
		//////////////
		//-- continue subdivision

		BvhTraversalData tFrontData;
		BvhTraversalData tBackData;
		
		// create new interior node and two leaf node
		currentNode = SubdivideNode(*sc, tFrontData, tBackData);

		// decrease the weighted average cost of the subdivisoin
		mTotalCost -= sc->GetRenderCostDecrease();
		mPvsEntries += sc->GetPvsEntriesIncr();

		// subdivision statistics
		if (1) PrintSubdivisionStats(*sc);


		///////////////////////////
		//-- push the new split candidates on the queue
		
		BvhSubdivisionCandidate *frontCandidate = 
				new BvhSubdivisionCandidate(tFrontData);
		BvhSubdivisionCandidate *backCandidate = 
				new BvhSubdivisionCandidate(tBackData);
		
		// preprocess view cells
		AssociateViewCellsWithObjects(*tData.mSampledObjects);

		EvalSubdivisionCandidate(*frontCandidate, true, false);
		EvalSubdivisionCandidate(*backCandidate, true, false);

		CollectDirtyCandidates(sc, dirtyList, true);
		ReleaseViewCells(*tData.mSampledObjects);
		
		// cross reference
		tFrontData.mNode->SetSubdivisionCandidate(frontCandidate); 
		tBackData.mNode->SetSubdivisionCandidate(backCandidate);

		//cout << "f: " << frontCandidate->GetPriority() << " b: " << backCandidate->GetPriority() << endl;
		tQueue.Push(frontCandidate);
		tQueue.Push(backCandidate);
	}

	/////////////////////////////////
	//-- node is a leaf => terminate traversal

	if (currentNode->IsLeaf())
	{
		/////////////////////
		//-- store additional info
		EvaluateLeafStats(tData);
        
		// this leaf is no candidate for splitting anymore
		// => detach subdivision candidate
		tData.mNode->SetSubdivisionCandidate(NULL); 
		// detach node so we don't delete it with the traversal data
		tData.mNode = NULL;
	}
	
	mSubdivTimer.Exit();

	return currentNode;
}


float BvHierarchy::EvalPriority(const BvhSubdivisionCandidate &splitCandidate,
								const float renderCostDecr,
								const float oldRenderCost) const
{
	float priority;

	if (mIsInitialSubdivision)
	{
		priority = (float)-splitCandidate.mParentData.mDepth;
		return priority;
	}

	BvhLeaf *leaf = splitCandidate.mParentData.mNode;

	// surface area heuristics is used when there is 
	// no view space subdivision available. 
	// In order to have some prioritized traversal, 
	// we use this formula instead
	if (mHierarchyManager->GetViewSpaceSubdivisionType() == 
		HierarchyManager::NO_VIEWSPACE_SUBDIV)
	{
		priority = EvalSahCost(leaf);
	}
	else
	{
		// take render cost of node into account 
		// otherwise danger of being stuck in a local minimum!
		const float factor = mRenderCostDecreaseWeight;

		priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;
		
		if (mHierarchyManager->mConsiderMemory)
		{
			priority /= ((float)splitCandidate.GetPvsEntriesIncr() + mMemoryConst);
		}
	}

	// hack: don't allow empty splits to be taken
	if (splitCandidate.mFrontObjects.empty() || splitCandidate.mBackObjects.empty())
		priority = 0;

	return priority;
}


static float AvgRayContribution(const int pvs, const int nRays)
{
	return (float)pvs / ((float)nRays + Limits::Small);
}


void BvHierarchy::EvalSubdivisionCandidate(BvhSubdivisionCandidate &splitCandidate, 
										   const bool computeSplitPlane,
										   const bool preprocessViewCells)
{
	mPlaneTimer.Entry();

#if STORE_VIEWCELLS_WITH_BVH
	if (preprocessViewCells) // fill view cells cache
		AssociateViewCellsWithObjects(*splitCandidate.mParentData.mSampledObjects);
#endif

	if (computeSplitPlane)
	{
		splitCandidate.mFrontObjects.clear();
		splitCandidate.mBackObjects.clear();
		splitCandidate.mSampledFrontObjects.clear();
		splitCandidate.mSampledBackObjects.clear();

		const bool sufficientSamples = 
			splitCandidate.mParentData.mNumRays > mMinRaysForVisibility;

		//if (!sufficientSamples) cout << splitCandidate.mParentData.mNumRays << " ";

		const bool useVisibiliyBasedHeuristics = 
					!mUseSah &&
					(mHierarchyManager->GetViewSpaceSubdivisionType() == 
					HierarchyManager::KD_BASED_VIEWSPACE_SUBDIV) &&
					sufficientSamples;

		// compute best object partition
		const float ratio =	SelectObjectPartition(splitCandidate.mParentData, 
												  splitCandidate.mFrontObjects, 
												  splitCandidate.mBackObjects,
												  useVisibiliyBasedHeuristics);
	
		// cost ratio violated?
		const bool maxCostRatioViolated = mTermMaxCostRatio < ratio;
		const int previousMisses = splitCandidate.mParentData.mMaxCostMisses;

		splitCandidate.SetMaxCostMisses(maxCostRatioViolated ? 
										previousMisses + 1 : previousMisses);

		StoreSampledObjects(splitCandidate.mSampledFrontObjects, splitCandidate.mFrontObjects);
		StoreSampledObjects(splitCandidate.mSampledBackObjects, splitCandidate.mBackObjects);
	}

	mPlaneTimer.Exit();

	mEvalTimer.Entry();

	const BvhTraversalData &tData = splitCandidate.mParentData;
	BvhLeaf *leaf = tData.mNode;

	// avg contribution of a ray to a pvs
	const float pvs = (float)CountViewCells(*tData.mSampledObjects);
	const float avgRayContri = AvgRayContribution((int)pvs, tData.mNumRays);

	// high avg ray contri, the result is influenced by undersampling
	splitCandidate.SetAvgRayContribution(avgRayContri);
	const float viewSpaceVol = GetViewSpaceVolume();

	const float oldVolume = EvalViewCellsVolume(*tData.mSampledObjects) / viewSpaceVol;
	const float oldRatio = (tData.mVolume) > 0 ? oldVolume / tData.mVolume : 1;
	const float parentVol = tData.mCorrectedVolume * oldRatio;

	// this leaf is a pvs entry in all the view cells
	// that see one of the objects.
	splitCandidate.mVolumeFrontViewCells = EvalViewCellsVolume(splitCandidate.mSampledFrontObjects) / viewSpaceVol;
	splitCandidate.mVolumeBackViewCells = EvalViewCellsVolume(splitCandidate.mSampledBackObjects) / viewSpaceVol;

	splitCandidate.mNumFrontViewCells = CountViewCells(splitCandidate.mSampledFrontObjects);
	splitCandidate.mNumBackViewCells = CountViewCells(splitCandidate.mSampledBackObjects);

	splitCandidate.mCorrectedFrontVolume = 
		mHierarchyManager->EvalCorrectedPvs(splitCandidate.mVolumeFrontViewCells, parentVol, avgRayContri);
	
	splitCandidate.mCorrectedBackVolume = 
		mHierarchyManager->EvalCorrectedPvs(splitCandidate.mVolumeBackViewCells, parentVol, avgRayContri);
	
	const float relfrontCost = splitCandidate.mCorrectedFrontVolume * 
		EvalAbsCost(splitCandidate.mFrontObjects);
	const float relBackCost =  splitCandidate.mCorrectedBackVolume * 
		EvalAbsCost(splitCandidate.mBackObjects);
	const float relParentCost = parentVol * EvalAbsCost(leaf->mObjects);

	// compute global decrease in render cost
	const float newRenderCost = relfrontCost + relBackCost;
	const float renderCostDecr = relParentCost - newRenderCost;
	
	splitCandidate.SetRenderCostDecrease(renderCostDecr);

	// increase in pvs entries
	const int pvsEntriesIncr = EvalPvsEntriesIncr(splitCandidate, 
												  avgRayContri, 
												  (int)pvs, 
												  splitCandidate.mNumFrontViewCells, 
												  splitCandidate.mNumBackViewCells);

	splitCandidate.SetPvsEntriesIncr(pvsEntriesIncr);

	if (0)
	{
		cout << "bvh volume cost"
		 << " avg ray contri: " << avgRayContri << " ratio: " << oldRatio 
		 << " parent: " << parentVol << " old vol: " << oldVolume 
		 << " frontvol: " << splitCandidate.mVolumeFrontViewCells << " corr. " << splitCandidate.mCorrectedFrontVolume 
		 << " backvol: " << splitCandidate.mVolumeBackViewCells << " corr. " << splitCandidate.mCorrectedBackVolume << endl;
	}

#ifdef GTP_DEBUG
	Debug << "old render cost: " << oldRenderCost << endl;
	Debug << "new render cost: " << newRenderCost << endl;
	Debug << "render cost decrease: " << renderCostDecr << endl;
#endif

	float priority = EvalPriority(splitCandidate, 
								  renderCostDecr,
								  relParentCost);

  	// compute global decrease in render cost
	splitCandidate.SetPriority(priority);

#if STORE_VIEWCELLS_WITH_BVH
	if (preprocessViewCells)
		ReleaseViewCells(*splitCandidate.mParentData.mSampledObjects);
#endif

	mEvalTimer.Exit();
}


int BvHierarchy::EvalPvsEntriesIncr(BvhSubdivisionCandidate &splitCandidate, 
									const float avgRayContri,
									const int numParentViewCells,
									const int numFrontViewCells,
									const int numBackViewCells) //const
{
	const float oldPvsSize = (float)numParentViewCells;//CountViewCells(*splitCandidate.mParentData.mSampledObjects);
	const float oldPvsRatio = (splitCandidate.mParentData.mPvs > 0) ? oldPvsSize / splitCandidate.mParentData.mPvs : 1;

	const float parentPvs = splitCandidate.mParentData.mCorrectedPvs * oldPvsRatio;

	const int frontViewCells = numFrontViewCells;//CountViewCells(splitCandidate.mSampledFrontObjects);
	const int backViewCells = numBackViewCells;//CountViewCells(splitCandidate.mSampledBackObjects);
	
	splitCandidate.mCorrectedFrontPvs = 
		mHierarchyManager->EvalCorrectedPvs((float)frontViewCells, parentPvs, avgRayContri);
	splitCandidate.mCorrectedBackPvs =
		mHierarchyManager->EvalCorrectedPvs((float)backViewCells, parentPvs, avgRayContri);

	if (0)
	cout << "bvh pvs"
		 << " avg ray contri: " << avgRayContri << " ratio: " << oldPvsRatio 
		 << " parent: " << parentPvs << " " << " old vol: " << oldPvsSize
		 << " frontpvs: " << frontViewCells << " corr. " << splitCandidate.mCorrectedFrontPvs 
		 << " backpvs: " << frontViewCells << " corr. " << splitCandidate.mCorrectedBackPvs << endl;

	return (int)(splitCandidate.mCorrectedFrontPvs + splitCandidate.mCorrectedBackPvs - parentPvs);
}


inline bool BvHierarchy::LocalTerminationCriteriaMet(const BvhTraversalData &tData) const
{
	const bool terminationCriteriaMet =
			(0
			|| ((int)tData.mNode->mObjects.size() <= 1)//mTermMinObjects)
			//|| (data.mProbability <= mTermMinProbability)
			//|| (data.mNumRays <= mTermMinRays)
		 );

#ifdef _DEBUG
	if (terminationCriteriaMet)
	{
		cout << "bvh local termination criteria met:" << endl;
		cout << "objects: " << (int)tData.mNode->mObjects.size() << " (" << mTermMinObjects << ")" << endl;
	}
#endif
	return terminationCriteriaMet; 
}


inline bool BvHierarchy::GlobalTerminationCriteriaMet(const BvhTraversalData &data) const
{
	// note: tracking for global cost termination 
	// does not make much sense for interleaved vsp / osp partition
	// as it is the responsibility of the hierarchy manager

	const bool terminationCriteriaMet =
		(0
		|| (mBvhStats.Leaves() >= mTermMaxLeaves)
		//|| (mBvhStats.mGlobalCostMisses >= mTermGlobalCostMissTolerance)
		//|| mOutOfMemory 
		);

#ifdef GTP_DEBUG
	if (terminationCriteriaMet)
	{
		Debug << "bvh global termination criteria met:" << endl;
		Debug << "cost misses: " << mBvhStats.mGlobalCostMisses << " " << mTermGlobalCostMissTolerance << endl;
		Debug << "leaves: " << mBvhStats.Leaves() << " " << mTermMaxLeaves << endl;
	}
#endif
	return terminationCriteriaMet; 
}


void BvHierarchy::EvaluateLeafStats(const BvhTraversalData &data)
{
	// the node became a leaf -> evaluate stats for leafs
	BvhLeaf *leaf = data.mNode;
	
	++ mCreatedLeaves;

	
	/*if (data.mProbability <= mTermMinProbability)
	{
		++ mBvhStats.minProbabilityNodes;
	}*/

	////////////////////////////////////////////
	// depth related stuff

	if (data.mDepth < mBvhStats.minDepth)
	{
		mBvhStats.minDepth = data.mDepth;
	}

	if (data.mDepth >= mTermMaxDepth)
	{
        ++ mBvhStats.maxDepthNodes;
	}

	// accumulate depth to compute average depth
	mBvhStats.accumDepth += data.mDepth;


	////////////////////////////////////////////
	// objects related stuff

	// note: the sum should alwaysbe total number of objects for bvh
	mBvhStats.objectRefs += (int)leaf->mObjects.size();

	if ((int)leaf->mObjects.size() <= mTermMinObjects)
	{
	     ++ mBvhStats.minObjectsNodes;
	}

	if (leaf->mObjects.empty())
	{
		++ mBvhStats.emptyNodes;
	}

	if ((int)leaf->mObjects.size() > mBvhStats.maxObjectRefs)
	{
		mBvhStats.maxObjectRefs = (int)leaf->mObjects.size();
	}

	if ((int)leaf->mObjects.size() < mBvhStats.minObjectRefs)
	{
		mBvhStats.minObjectRefs = (int)leaf->mObjects.size();
	}

	////////////////////////////////////////////
	// ray related stuff

	// note: this number should always accumulate to the total number of rays
	mBvhStats.rayRefs += data.mNumRays;
	
	if (data.mNumRays <= mTermMinRays)
	{
	     ++ mBvhStats.minRaysNodes;
	}

	if (data.mNumRays > mBvhStats.maxRayRefs)
	{
		mBvhStats.maxRayRefs = data.mNumRays;
	}

	if (data.mNumRays < mBvhStats.minRayRefs)
	{
		mBvhStats.minRayRefs = data.mNumRays;
	}

#ifdef _DEBUG
	cout << "depth: " << data.mDepth << " objects: " << (int)leaf->mObjects.size() 
		 << " rays: " << data.mNumRays << " rays / objects " 
		 << (float)data.mNumRays / (float)leaf->mObjects.size() << endl;
#endif
}


#if 1

/// compute object boundaries using spatial mid split
float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
											const int axis,
											ObjectContainer &objectsFront,
											ObjectContainer &objectsBack)
{
	AxisAlignedBox3 parentBox = tData.mNode->GetBoundingBox();

	const float maxBox = parentBox.Max(axis);
	const float minBox = parentBox.Min(axis);

	float midPoint = (maxBox + minBox) * 0.5f;

	ObjectContainer::const_iterator oit, oit_end = tData.mNode->mObjects.end();
	
	for (oit = tData.mNode->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		const AxisAlignedBox3 box = obj->GetBox();

		const float objMid = (box.Max(axis) + box.Min(axis)) * 0.5f;

		// object mailed => belongs to back objects
		if (objMid < midPoint) 
		{
			objectsBack.push_back(obj);
		}
		else
		{
			objectsFront.push_back(obj);
		}
	}

	AxisAlignedBox3 fbox = EvalBoundingBox(objectsFront, &parentBox);
	AxisAlignedBox3 bbox = EvalBoundingBox(objectsBack, &parentBox);

	const float oldRenderCost = (float)tData.mNode->mObjects.size() * parentBox.SurfaceArea();
	const float newRenderCost = (float)objectsFront.size() * fbox.SurfaceArea() +  (float)objectsBack.size() * bbox.SurfaceArea();

	const float ratio = newRenderCost / oldRenderCost;
	return ratio;
}

#else

/// compute object partition by getting balanced objects on the left and right side
float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
											const int axis,
											ObjectContainer &objectsFront,
											ObjectContainer &objectsBack)
{
	PrepareLocalSubdivisionCandidates(tData, axis);
	
	SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();

	int i = 0;
	const int border = (int)tData.mNode->mObjects.size() / 2;

    for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ i)
	{
		Intersectable *obj = (*cit).mObject;

		// object mailed => belongs to back objects
		if (i < border) 
		{
			objectsBack.push_back(obj);
		}
		else
		{
			objectsFront.push_back(obj);
		}
	}

#if 1
	// hack: always take driving axis
	const float cost = (tData.mNode->GetBoundingBox().Size().DrivingAxis() == axis) ? -1.0f : 0.0f;
#else
	const float oldRenderCost = EvalAbsCost(tData.mLeaf->mObjects) / EvalProbability(tData.mSampledObjects);
	const float newRenderCost = EvalRenderCost(objectsFront) + EvalRenderCost(objectsBack);

	const float cost = newRenderCost / oldRenderCost;
#endif

	return cost;
}
#endif

#if 1

float BvHierarchy::EvalSah(const BvhTraversalData &tData,
						   const int axis,
						   ObjectContainer &objectsFront,
						   ObjectContainer &objectsBack)
{
	// go through the lists, count the number of objects left and right
	// and evaluate the following cost funcion:
	// C = ct_div_ci  + (ol + or) / queries
	PrepareLocalSubdivisionCandidates(tData, axis);

	const float totalRenderCost = EvalAbsCost(tData.mNode->mObjects);
	float objectsLeft = 0, objectsRight = totalRenderCost;
  
	const AxisAlignedBox3 nodeBbox = tData.mNode->GetBoundingBox();
	const float boxArea = nodeBbox.SurfaceArea();

	float minSum = 1e20f;
  
	float minBorder = nodeBbox.Max(axis);
	float maxBorder = nodeBbox.Min(axis);

	float areaLeft = 0, areaRight = 0;

	SortableEntryContainer::const_iterator currentPos =
		mSubdivisionCandidates->begin();
	
	vector<float> bordersRight;

	// we keep track of both borders of the bounding boxes =>
	// store the events in descending order

	bordersRight.resize(mSubdivisionCandidates->size());

	SortableEntryContainer::reverse_iterator rcit = 
		mSubdivisionCandidates->rbegin(), rcit_end = mSubdivisionCandidates->rend();

	vector<float>::reverse_iterator rbit = bordersRight.rbegin();

	for (; rcit != rcit_end; ++ rcit, ++ rbit) 
	{
		Intersectable *obj = (*rcit).mObject;
		const AxisAlignedBox3 obox = obj->GetBox();

		if (obox.Min(axis) < minBorder)
		{
			minBorder = obox.Min(axis);
		}

		(*rbit) = minBorder;
	}

	// record surface areas during the sweep
	float al = 0;
	float ar = boxArea;

	vector<float>::const_iterator bit = bordersRight.begin();
	SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();

	for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit) 
	{
		Intersectable *obj = (*cit).mObject;

#if USE_BETTER_RENDERCOST_EST
		const float renderCost = ViewCellsManager::EvalRenderCost(obj);
		
		objectsLeft += renderCost;
		objectsRight -= renderCost;

		const bool noValidSplit = ((objectsLeft <= Limits::Small) || (objectsRight <= Limits::Small));

#else
		++ objectsLeft;
		-- objectsRight;

		const bool noValidSplit = !objectsLeft || !objectsRight;
#endif
		const AxisAlignedBox3 obox = obj->GetBox();

		// the borders of the bounding boxes have changed
		if (obox.Max(axis) > maxBorder)
		{
			maxBorder = obox.Max(axis);
		}

		minBorder = (*bit);

		AxisAlignedBox3 lbox = nodeBbox;
		AxisAlignedBox3 rbox = nodeBbox;

		lbox.SetMax(axis, maxBorder);
		rbox.SetMin(axis, minBorder);

		al = lbox.SurfaceArea();
		ar = rbox.SurfaceArea();

		const float sum = noValidSplit ? 1e25f : objectsLeft * al + objectsRight * ar;
      
		/*cout << "pos=" << (*cit).mPos << "\t q=(" << objectsLeft << "," << objectsRight <<")\t r=(" 
			 << lbox.SurfaceArea() << "," << rbox.SurfaceArea() << ")" << endl;
		cout << "minborder: " << minBorder << " maxborder: " << maxBorder << endl;
	    cout << "cost= " << sum << endl;*/
	
		if (sum < minSum) 
		{	
			minSum = sum;
			areaLeft = al;
			areaRight = ar;

			// objects belong to left side now
			for (; currentPos != (cit + 1); ++ currentPos);
		}
	}

	////////////
	//-- assign object to front and back volume

	// belongs to back bv
	for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
		objectsBack.push_back((*cit).mObject);
	
	// belongs to front bv
	for (cit = currentPos; cit != cit_end; ++ cit)
		objectsFront.push_back((*cit).mObject);

	float newCost = minSum / boxArea;
	float ratio = newCost / totalRenderCost;
  
#ifdef GTP_DEBUG
	cout << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of " 
		 << (int)tData.mNode->mObjects.size() << ")\t area=(" 
		 << areaLeft << ", " << areaRight << ", " << boxArea << ")" << endl
		 << "cost= " << newCost << " oldCost=" << totalRenderCost / boxArea << endl;
#endif

	return ratio;
}

#else

float BvHierarchy::EvalSah(const BvhTraversalData &tData,
						   const int axis,
						   ObjectContainer &objectsFront,
						   ObjectContainer &objectsBack)
{
	// go through the lists, count the number of objects left and right
	// and evaluate the following cost funcion:
	// C = ct_div_ci  + (ol + or) / queries
	PrepareLocalSubdivisionCandidates(tData, axis);

	const float totalRenderCost = EvalAbsCost(tData.mNode->mObjects);
	float objectsLeft = 0, objectsRight = totalRenderCost;
  
	const AxisAlignedBox3 nodeBbox = tData.mNode->GetBoundingBox();

	const float minBox = nodeBbox.Min(axis);
	const float maxBox = nodeBbox.Max(axis);
	const float boxArea = nodeBbox.SurfaceArea();

	float minSum = 1e20f;
  
	Vector3 minBorder = nodeBbox.Max();
	Vector3 maxBorder = nodeBbox.Min();

	float areaLeft = 0, areaRight = 0;

	SortableEntryContainer::const_iterator currentPos =
		mSubdivisionCandidates->begin();
	
	vector<Vector3> bordersRight;

	// we keep track of both borders of the bounding boxes =>
	// store the events in descending order
	bordersRight.resize(mSubdivisionCandidates->size());

	SortableEntryContainer::reverse_iterator rcit = 
		mSubdivisionCandidates->rbegin(), rcit_end = 
		mSubdivisionCandidates->rend();

	vector<Vector3>::reverse_iterator rbit = bordersRight.rbegin();

	for (; rcit != rcit_end; ++ rcit, ++ rbit) 
	{
		Intersectable *obj = (*rcit).mObject;
		const AxisAlignedBox3 obox = obj->GetBox();

		for (int i = 0; i < 3; ++ i)
		{
			if (obox.Min(i) < minBorder[i])
			{
				minBorder[i] = obox.Min(i);
			}
		}

		(*rbit) = minBorder;
	}

	// temporary surface areas
	float al = 0;
	float ar = boxArea;

	vector<Vector3>::const_iterator bit = bordersRight.begin();
	SortableEntryContainer::const_iterator cit, cit_end =
		mSubdivisionCandidates->end();

	for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit) 
	{
		Intersectable *obj = (*cit).mObject;

		const float renderCost = ViewCellsManager::EvalRenderCost(obj);
		
		objectsLeft += renderCost;
		objectsRight -= renderCost;

		const AxisAlignedBox3 obox = obj->GetBox();

		AxisAlignedBox3 lbox = nodeBbox;
		AxisAlignedBox3 rbox = nodeBbox;
	
		// the borders of the left bounding box have changed
		for (int i = 0; i < 3; ++ i)
		{
			if (obox.Max(i) > maxBorder[i])
			{
				maxBorder[i] = obox.Max(i);
			}
		}

		minBorder = (*bit);

		lbox.SetMax(maxBorder);
		rbox.SetMin(minBorder);

		al = lbox.SurfaceArea();
		ar = rbox.SurfaceArea();
	
		const bool noValidSplit = ((objectsLeft <= Limits::Small) || (objectsRight <= Limits::Small));
		const float sum =  noValidSplit ? 1e25 : objectsLeft * al + objectsRight * ar;
      
		/*cout << "pos=" << (*cit).mPos << "\t q=(" << objectsLeft << "," << objectsRight <<")\t r=(" 
			 << lbox.SurfaceArea() << "," << rbox.SurfaceArea() << ")" << endl;
		cout << "minborder: " << minBorder << " maxborder: " << maxBorder << endl;
	    cout << "cost= " << sum << endl;*/
	
		if (sum < minSum) 
		{	
			minSum = sum;
			areaLeft = al;
			areaRight = ar;

			// objects belong to left side now
			for (; currentPos != (cit + 1); ++ currentPos);
		}
	}

	/////////////
	//-- assign object to front and back volume

	// belongs to back bv
	for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
		objectsBack.push_back((*cit).mObject);
	
	// belongs to front bv
	for (cit = currentPos; cit != cit_end; ++ cit)
		objectsFront.push_back((*cit).mObject);

	float newCost = minSum / boxArea;
	float ratio = newCost / totalRenderCost;
  
#ifdef GTP_DEBUG
	cout << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of " 
		 << (int)tData.mNode->mObjects.size() << ")\t area=(" 
		 << areaLeft << ", " << areaRight << ", " << boxArea << ")" << endl
		 << "cost= " << newCost << " oldCost=" << totalRenderCost / boxArea << endl;
#endif

	return ratio;
}

#endif

static bool PrepareOutput(const int axis,
						  const int leaves,
						  ofstream &sumStats, 
						  ofstream &vollStats, 
						  ofstream &volrStats)
{
	if ((axis == 0) && (leaves > 0) && (leaves < 90))
	{
		char str[64]; 	
		sprintf(str, "tmp/bvh_heur_sum-%04d.log", leaves);
		sumStats.open(str);
		sprintf(str, "tmp/bvh_heur_voll-%04d.log", leaves);
		vollStats.open(str);
		sprintf(str, "tmp/bvh_heur_volr-%04d.log", leaves);
		volrStats.open(str);
	}

	return sumStats.is_open() && vollStats.is_open() && volrStats.is_open();
}


static void PrintHeuristics(const float objectsRight,
							const float sum,
							const float volLeft,
							const float volRight,
							const float viewSpaceVol,
							ofstream &sumStats, 
							ofstream &vollStats, 
							ofstream &volrStats)
{
	sumStats
		<< "#Position\n" << objectsRight << endl
		<< "#Sum\n" << sum / viewSpaceVol << endl
		<< "#Vol\n" << (volLeft +  volRight) / viewSpaceVol << endl;

	vollStats
		<< "#Position\n" << objectsRight << endl
		<< "#Vol\n" << volLeft / viewSpaceVol << endl;

	volrStats
		<< "#Position\n" << objectsRight << endl
		<< "#Vol\n" << volRight / viewSpaceVol << endl;
}


float BvHierarchy::EvalLocalCostHeuristics(const BvhTraversalData &tData,
										   const int axis,
										   ObjectContainer &objectsFront,
										   ObjectContainer &objectsBack)
{
	/////////////////////////////////////////////
	//-- go through the lists, count the number of objects 
	//-- left and right and evaluate the cost funcion

	// prepare the heuristics by setting mailboxes and counters
	const float totalVol = PrepareHeuristics(tData, axis);
	
	// local helper variables
	float volLeft = 0;
	float volRight = totalVol;
	
#if 1//USE_BETTER_RENDERCOST_EST
	const float nTotalObjects = EvalAbsCost(tData.mNode->mObjects);
	float nObjectsLeft = 0;
	float nObjectsRight = nTotalObjects;
#else
	const int nTotalObjects = (int)EvalAbsCost(tData.mNode->mObjects);
	int nObjectsLeft = 0;
	int nObjectsRight = (int)nTotalObjects;
#endif

	const float viewSpaceVol = 
		mViewCellsManager->GetViewSpaceBox().GetVolume();

	SortableEntryContainer::const_iterator backObjectsStart = 
		mSubdivisionCandidates->begin();

	/////////////////////////////////
	//-- the parameters for the current optimum

	float volBack = volLeft;
	float volFront = volRight;
	float newRenderCost = nTotalObjects * totalVol;

#ifdef GTP_DEBUG
	ofstream sumStats;
	ofstream vollStats;
	ofstream volrStats;

	const bool printStats = PrepareOutput(axis, 
										  mBvhStats.Leaves(), 
										  sumStats, 
										  vollStats, 
										  volrStats);
#endif

	///////////////////////
	//-- the sweep heuristics
	//-- traverse through events and find best split plane

	SortableEntryContainer::const_iterator cit, 
		cit_end = cit_end = mSubdivisionCandidates->end();

	for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit)
	{
		Intersectable *object = (*cit).mObject;
	
		// evaluate change in l and r volume
		// voll = view cells that see only left node (i.e., left pvs)
		// volr = view cells that see only right node (i.e., right pvs)
		EvalHeuristicsContribution(object, volLeft, volRight);

#if USE_BETTER_RENDERCOST_EST

		const float rc = ViewCellsManager::EvalRenderCost(object);
		
		nObjectsLeft += rc;
		nObjectsRight -= rc;
	
#else
		++ nObjectsLeft;
		-- nObjectsRight;
#endif
	
		// split is only valid if #objects on left and right is not zero
		const bool noValidSplit = (nObjectsRight <= Limits::Small);

		// the heuristics
	    const float sum = noValidSplit ? 
			1e25f : volLeft * (float)nObjectsLeft + volRight * (float)nObjectsRight;

#ifdef GTP_DEBUG
		if (printStats)
		{
			PrintHeuristics(nObjectsRight, sum, volLeft, 
							volRight, viewSpaceVol,
							sumStats, vollStats, volrStats);
		}
#endif

		if (sum < newRenderCost)
		{
			newRenderCost = sum;

			volBack = volLeft;
			volFront = volRight;

			// objects belongs to left side now
			for (; backObjectsStart != (cit + 1); ++ backObjectsStart);
		}
	}

	////////////////////////////////////////
	//-- assign object to front and back volume

	// belongs to back bv
	for (cit = mSubdivisionCandidates->begin(); cit != backObjectsStart; ++ cit)
	{
		objectsBack.push_back((*cit).mObject);
	}
	// belongs to front bv
	for (cit = backObjectsStart; cit != cit_end; ++ cit)
	{
		objectsFront.push_back((*cit).mObject);
	}

	// render cost of the old parent
	const float oldRenderCost = (float)nTotalObjects * totalVol + Limits::Small;
	// the relative cost ratio
	const float ratio = newRenderCost / oldRenderCost;

#ifdef GTP_DEBUG
	Debug << "\n§§§§ bvh eval const decrease §§§§" << endl
		  << "back pvs: " << (int)objectsBack.size() << " front pvs: " 
		  << (int)objectsFront.size() << " total pvs: " << nTotalObjects << endl 
		  << "back p: " << volBack / viewSpaceVol << " front p " 
		  << volFront / viewSpaceVol << " p: " << totalVol / viewSpaceVol << endl
		  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " 
		  << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " 
		  << oldRenderCost / viewSpaceVol - newRenderCost / viewSpaceVol << endl;
#endif

	return ratio;
}


void BvHierarchy::PrepareLocalSubdivisionCandidates(const BvhTraversalData &tData,
													const int axis)											
{
	mSortTimer.Entry();
	
	//-- insert object queries
	ObjectContainer *objects = mUseGlobalSorting ? 
		tData.mSortedObjects[axis] : &tData.mNode->mObjects;

	CreateLocalSubdivisionCandidates(*objects, &mSubdivisionCandidates, !mUseGlobalSorting, axis);
	
	mSortTimer.Exit();
}


void BvHierarchy::CreateLocalSubdivisionCandidates(const ObjectContainer &objects, 
												  SortableEntryContainer **subdivisionCandidates, 
												  const bool sortEntries,
												  const int axis)
{
	(*subdivisionCandidates)->clear();

	// compute requested size and look if subdivision candidate has to be recomputed
	const int requestedSize = (int)objects.size();
	
	// creates a sorted split candidates array
	if ((*subdivisionCandidates)->capacity() > 500000 &&
		requestedSize < (int)((*subdivisionCandidates)->capacity() / 10) )
	{
    	delete (*subdivisionCandidates);
   		(*subdivisionCandidates) = new SortableEntryContainer;
	}

	(*subdivisionCandidates)->reserve(requestedSize);

	ObjectContainer::const_iterator oit, oit_end = objects.end();

	for (oit = objects.begin(); oit < oit_end; ++ oit)
	{
		(*subdivisionCandidates)->push_back(SortableEntry(*oit, (*oit)->GetBox().Center(axis)));
	}

	if (sortEntries)
	{	// no presorted candidate list
		stable_sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
		//sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
	}
}


const BvhStatistics &BvHierarchy::GetStatistics() const
{
	return mBvhStats;
}


float BvHierarchy::PrepareHeuristics(const BvhTraversalData &tData, 
									 const int axis)
{	
	BvhLeaf *leaf = tData.mNode;
	float vol = 0;

    // sort so we can use a sweep from right to left
	PrepareLocalSubdivisionCandidates(tData, axis);
	
	// collect and mark the view cells as belonging to front pvs
	ViewCellContainer viewCells;

	const bool setCounter = true;
	const bool onlyUnmailed = true;

	
	CollectViewCells(*tData.mSampledObjects, 
					 viewCells, 
					 setCounter, 
					 onlyUnmailed);

	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
#if USE_VOLUMES_FOR_HEURISTICS
		const float volIncr = (*vit)->GetVolume();
#else
		const float volIncr = 1.0f;
#endif
		vol += volIncr;
	}

	// mail view cells that go from front node to back node
	ViewCell::NewMail();
	
	return vol;
}

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


void BvHierarchy::EvalHeuristicsContribution(Intersectable *obj,
											 float &volLeft, 
											 float &volRight)
{
	// collect all view cells associated with this objects 
	// (also multiple times, if they are pierced by several rays)
	ViewCellContainer viewCells;

	const bool useMailboxing = false;
	const bool setCounter = false;
	const bool onlyUnmailedRays = true;

	CollectViewCells(obj, viewCells, useMailboxing, setCounter, onlyUnmailedRays);

	// classify view cells and compute volume contri accordingly
	// possible view cell classifications:
	// view cell mailed => view cell can be seen from left child node
	// view cell counter > 0 view cell can be seen from right child node
	// combined: view cell volume belongs to both nodes
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
	
	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
		// view cells can also be seen from left child node
		ViewCell *viewCell = *vit;

#if USE_VOLUMES_FOR_HEURISTICS
		const float vol = viewCell->GetVolume();
#else
		const float vol = 1.0f;
#endif
		if (!viewCell->Mailed())
		{
			viewCell->Mail();
			// we now see view cell from both nodes 
			// => add volume to left node
			volLeft += vol;
		}

		// last reference into the right node
		if (-- viewCell->mCounter == 0)
		{	
			// view cell was previously seen from both nodes  =>
			// remove volume from right node
			volRight -= vol;
		}
	}
}


void BvHierarchy::SetViewCellsManager(ViewCellsManager *vcm)
{
	mViewCellsManager = vcm;
}


AxisAlignedBox3 BvHierarchy::GetBoundingBox() const
{
	return mBoundingBox;
}


float BvHierarchy::SelectObjectPartition(const BvhTraversalData &tData,
										 ObjectContainer &frontObjects,
										 ObjectContainer &backObjects,
										 bool useVisibilityBasedHeuristics)
{
	mSplitTimer.Entry();

	if (mIsInitialSubdivision)
	{
		ApplyInitialSplit(tData, frontObjects, backObjects);
		return 0;
	}

	ObjectContainer nFrontObjects[3];
	ObjectContainer nBackObjects[3];
	float nCostRatio[3];

	int sAxis = 0;
	int bestAxis = -1;

	if (mOnlyDrivingAxis)
	{
		const AxisAlignedBox3 box = tData.mNode->GetBoundingBox();
		sAxis = box.Size().DrivingAxis();
	}

	// if #rays high consider only use a subset of the rays for 
	// visibility based heuristics
	VssRay::NewMail();

	/*if ((mMaxTests < tData.mNumRays) && mUseCostHeuristics && useVisibilityBasedHeuristics)
	{
		VssRayContainer rays;

		// maximal 2 objects share the same ray
		rays.reserve(tData.mNumRays * 2);
		CollectRays(tData.mNode->mObjects, rays);

		const float prop = (float)mMaxTests / (float)rays.size();

		VssRayContainer::const_iterator rit, rit_end = rays.end();

		// mail rays which will not be considered
		for (rit = rays.begin(); rit != rit_end; ++ rit)
		{
			if (Random(1.0f) > prop)
			{
				(*rit)->Mail();
			}
		}		
	}*/

	////////////////////////////////////
	//-- evaluate split cost for all three axis
	
	for (int axis = 0; axis < 3; ++ axis)
	{
		if (!mOnlyDrivingAxis || (axis == sAxis))
		{
			if (mUseCostHeuristics)
			{
				//////////////////////////////////
                //-- split objects using heuristics
				
				if (useVisibilityBasedHeuristics)
				{
					///////////
					//-- heuristics using objects weighted by view cells volume
					nCostRatio[axis] =
						EvalLocalCostHeuristics(tData, 
												axis, 
												nFrontObjects[axis], 
												nBackObjects[axis]);
				}
				else
				{	
					//////////////////
					//-- view cells not constructed yet	=> use surface area heuristic			
					nCostRatio[axis] = EvalSah(tData, 
											   axis, 
											   nFrontObjects[axis], 
											   nBackObjects[axis]);
				}
			}
			else
			{
				//-- split objects using some simple criteria
				nCostRatio[axis] =
					EvalLocalObjectPartition(tData, axis, nFrontObjects[axis], nBackObjects[axis]);
			}

			// avoid splits in degenerate axis with high penalty
			if (1 &&
				(tData.mNode->GetBoundingBox().Size(axis) < 0.0001))//Limits::Small))
			{
				nCostRatio[axis] += 9999;
			}

			if ((bestAxis == -1) || (nCostRatio[axis] < nCostRatio[bestAxis]))
			{
				bestAxis = axis;
			}
		}
	}

    ////////////////
	//-- assign values

	frontObjects = nFrontObjects[bestAxis];
	backObjects = nBackObjects[bestAxis];

	mSplitTimer.Exit();

	//cout << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
	return nCostRatio[bestAxis];
}


int BvHierarchy::AssociateObjectsWithRays(const VssRayContainer &rays) const
{
	int nRays = 0;
	VssRayContainer::const_iterator rit, rit_end = rays.end();

	VssRay::NewMail();

    for (rit = rays.begin(); rit != rays.end(); ++ rit)
	{
		VssRay *ray = (*rit);

		if (ray->mTerminationObject)
		{
			ray->mTerminationObject->GetOrCreateRays()->push_back(ray);
			if (!ray->Mailed())
			{
				ray->Mail();
				++ nRays;
			}
		}

#if COUNT_ORIGIN_OBJECTS

		if (ray->mOriginObject)
		{
			ray->mOriginObject->GetOrCreateRays()->push_back(ray);

			if (!ray->Mailed())
			{
				ray->Mail();
				++ nRays;
			}
		}
#endif
	}

	return nRays;
}


void BvHierarchy::PrintSubdivisionStats(const SubdivisionCandidate &sc)
{
	const float costDecr = sc.GetRenderCostDecrease();	

	mSubdivisionStats 
			<< "#Leaves\n" << mBvhStats.Leaves() << endl
			<< "#RenderCostDecrease\n" << costDecr << endl 
			<< "#TotalRenderCost\n" << mTotalCost << endl
			<< "#EntriesInPvs\n" << mPvsEntries << endl;
}


void BvHierarchy::CollectRays(const ObjectContainer &objects,
							  VssRayContainer &rays) const
{
	VssRay::NewMail();
	ObjectContainer::const_iterator oit, oit_end = objects.end();

	// evaluate reverse pvs and view cell volume on left and right cell
	// note: should I take all leaf objects or rather the objects hit by rays?
	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

		for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
		{
			VssRay *ray = (*rit);

			if (!ray->Mailed())
			{
				ray->Mail();
				rays.push_back(ray);
			}
		}
	}
}


float BvHierarchy::EvalAbsCost(const ObjectContainer &objects)
{
#if USE_BETTER_RENDERCOST_EST
	ObjectContainer::const_iterator oit, oit_end = objects.end();

	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		objRenderCost += ViewCellsManager::GetRendercost(*oit);
	}
#else
	return (float)objects.size();
#endif
}


float BvHierarchy::EvalSahCost(BvhLeaf *leaf) const
{
	////////////////
	//-- surface area heuristics

	if (leaf->mObjects.empty())
		return 0.0f;

	const AxisAlignedBox3 box = GetBoundingBox(leaf);
	const float area = box.SurfaceArea();
	const float viewSpaceArea = mViewCellsManager->GetViewSpaceBox().SurfaceArea();

	return EvalAbsCost(leaf->mObjects) * area / viewSpaceArea;
}


float BvHierarchy::EvalRenderCost(const ObjectContainer &objects)// const
{	
	///////////////
	//-- render cost heuristics
	const float objRenderCost = EvalAbsCost(objects);

	const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();

	// probability that view point lies in a view cell which sees this node
	const float p = EvalViewCellsVolume(objects) / viewSpaceVol;
	
	return objRenderCost * p;
}


float BvHierarchy::EvalProbability(const ObjectContainer &objects)// const
{	
	const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
	
	// probability that view point lies in a view cell which sees this node
	return EvalViewCellsVolume(objects) / viewSpaceVol;
}


AxisAlignedBox3 BvHierarchy::EvalBoundingBox(const ObjectContainer &objects,
											 const AxisAlignedBox3 *parentBox) const
{
	// if there are no objects in this box, box size is set to parent box size. 
	// Question: Invalidate box instead?
	if (parentBox && objects.empty())
		return *parentBox;

	AxisAlignedBox3 box;
	box.Initialize();

	ObjectContainer::const_iterator oit, oit_end = objects.end();

	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		// grow bounding box to include all objects
		box.Include(obj->GetBox());
	}

	return box;
}


void BvHierarchy::CollectLeaves(BvhNode *root, vector<BvhLeaf *> &leaves) const
{
	stack<BvhNode *> nodeStack;
	nodeStack.push(root);

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

		if (node->IsLeaf()) 
		{
			BvhLeaf *leaf = (BvhLeaf *)node;
			leaves.push_back(leaf);
		}
		else 
		{
			BvhInterior *interior = (BvhInterior *)node;

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


void BvHierarchy::CollectNodes(BvhNode *root, vector<BvhNode *> &nodes) const
{
	stack<BvhNode *> nodeStack;
	nodeStack.push(root);

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

		nodes.push_back(node);
		
		if (!node->IsLeaf()) 
		{
			BvhInterior *interior = (BvhInterior *)node;

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


AxisAlignedBox3 BvHierarchy::GetBoundingBox(BvhNode *node) const
{
	return node->GetBoundingBox();
}


int BvHierarchy::CollectViewCells(const ObjectContainer &objects, 
								  ViewCellContainer &viewCells,
								  const bool setCounter,
								  const bool onlyUnmailedRays)// const
{
	ViewCell::NewMail();

	ObjectContainer::const_iterator oit, oit_end = objects.end();

	// use mailing to avoid dublicates
	const bool useMailBoxing = true;

	int numRays = 0;
	// loop through all object and collect view cell pvs of this node
	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		// use mailing to avoid duplicates
		numRays += CollectViewCells(*oit, viewCells, useMailBoxing, setCounter, onlyUnmailedRays);
	}

	return numRays;
}


#if STORE_VIEWCELLS_WITH_BVH


void BvHierarchy::ReleaseViewCells(const ObjectContainer &objects)
{
	ObjectContainer::const_iterator oit, oit_end = objects.end();

	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		(*oit)->DelViewCells();
	}
}


void BvHierarchy::AssociateViewCellsWithObjects(const ObjectContainer &objects) const
{
	ObjectContainer::const_iterator oit, oit_end = objects.end();

	const bool useMailBoxing = true;
	VssRay::NewMail();
	
	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
			ViewCell::NewMail();
			// use mailing to avoid duplicates
			AssociateViewCellsWithObject(*oit, useMailBoxing);
	}
}


int BvHierarchy::AssociateViewCellsWithObject(Intersectable *obj, const bool useMailBoxing) const
{
	int nRays = 0;

	if (!obj->GetOrCreateViewCells()->empty())
	{
		cerr << "AssociateViewCellsWithObject: view cells cache not working" << endl;
	}

	ViewCellContainer *objViewCells = obj->GetOrCreateViewCells();
	VssRayContainer *vssRays = obj->GetOrCreateRays();

	VssRayContainer::const_iterator rit, rit_end = vssRays->end();

	// fill cache
	for (rit = vssRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);

		//	if (onlyUnmailedRays && ray->Mailed())
		//		continue;
		mHierarchyManager->mVspTree->GetViewCells(*ray, *objViewCells);
		
		if (!useMailBoxing || !ray->Mailed())
		{
			if (useMailBoxing)
				ray->Mail();

			++ nRays;
		}
	}

	return nRays;
}



int BvHierarchy::CountViewCells(Intersectable *obj) //const
{
	ViewCellContainer *viewCells = obj->GetOrCreateViewCells();

	if (obj->GetOrCreateViewCells()->empty())
	{
		//cerr << "h";//CountViewCells: view cells empty, view cells cache not working" << endl;
		return CountViewCellsFromRays(obj);
	}
	
	int result = 0;

	ViewCellContainer::const_iterator vit, vit_end = viewCells->end();
	
	for (vit = viewCells->begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;

		// store view cells
		if (!vc->Mailed())
		{
			vc->Mail();
			++ result;
		}
	}

	return result;
}


int BvHierarchy::CollectViewCells(Intersectable *obj, 
								  ViewCellContainer &viewCells,
								  const bool useMailBoxing,
								  const bool setCounter,
								  const bool onlyUnmailedRays)// const
{
	if (obj->GetOrCreateViewCells()->empty())
	{
		//cerr << "g";//CollectViewCells: view cells empty, view cells cache not working" << endl;
		return CollectViewCellsFromRays(obj, viewCells, useMailBoxing, setCounter, onlyUnmailedRays);
	}

	mCollectTimer.Entry();

	ViewCellContainer *objViewCells = obj->GetOrCreateViewCells();

	///////////
	//-- use view cells cache

	// fastest way: just copy to the end
	//if (!useMailBoxing)	viewCells.insert(viewCells.end(), objViewCells->begin(), objViewCells->end());

	// loop through view cells
	// matt: probably slow to insert  view cells one by one
	ViewCellContainer::const_iterator vit, vit_end = objViewCells->end();

	for (vit = objViewCells->begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;

		// store view cells
		if (!useMailBoxing || !vc->Mailed())
		{
			if (useMailBoxing) 
			{
				// view cell not mailed
				vc->Mail();
				
				if (setCounter)
					vc->mCounter = 0;
				//viewCells.push_back(vc);
			}

			viewCells.push_back(vc);
		}

		if (setCounter)
			++ vc->mCounter;
	}

	mCollectTimer.Exit();

	return (int)objViewCells->size();
}


int BvHierarchy::CollectViewCellsFromRays(Intersectable *obj, 
										  ViewCellContainer &viewCells,
										  const bool useMailBoxing,
										  const bool setCounter,
										  const bool onlyUnmailedRays)
{
	mCollectTimer.Entry();
	VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

	int numRays = 0;

	for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);

		if (onlyUnmailedRays && ray->Mailed())
			continue;
		
		++ numRays;

		ViewCellContainer tmpViewCells;
		mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);

		// matt: probably slow to allocate memory for view cells every time
		ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();

		for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// store view cells
			if (!useMailBoxing || !vc->Mailed())
			{
				if (useMailBoxing) // => view cell not mailed
				{
					vc->Mail();
					if (setCounter)
						vc->mCounter = 0;
				}

				viewCells.push_back(vc);
			}
			
			if (setCounter)
				++ vc->mCounter;
		}
	}

	mCollectTimer.Exit();
	return numRays;
}


int BvHierarchy::CountViewCellsFromRays(Intersectable *obj) //const
{
	int result = 0;
	
	VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

	for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);
		ViewCellContainer tmpViewCells;
	
		mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
		
		ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
		for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// store view cells
			if (!vc->Mailed())
			{
				vc->Mail();
				++ result;
			}
		}
	}

	return result;
}

#else

int BvHierarchy::CountViewCells(Intersectable *obj) //const
{
	int result = 0;
	
	VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

	for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);
		ViewCellContainer tmpViewCells;
	
		mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
		
		ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
		for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// store view cells
			if (!vc->Mailed())
			{
				vc->Mail();
				++ result;
			}
		}
	}

	return result;
}


int BvHierarchy::CollectViewCells(Intersectable *obj, 
								  ViewCellContainer &viewCells,
								  const bool useMailBoxing,
								  const bool setCounter,
								  const bool onlyUnmailedRays)
{
	mCollectTimer.Entry();
	VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

	int numRays = 0;

	for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);

		if (onlyUnmailedRays && ray->Mailed())
			continue;
		
		++ numRays;

		ViewCellContainer tmpViewCells;
		mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);

		// matt: probably slow to allocate memory for view cells every time
		ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();

		for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// store view cells
			if (!useMailBoxing || !vc->Mailed())
			{
				if (useMailBoxing) // => view cell not mailed
				{
					vc->Mail();
					if (setCounter)
						vc->mCounter = 0;
				}

				viewCells.push_back(vc);
			}
			
			if (setCounter)
				++ vc->mCounter;
		}
	}

	mCollectTimer.Exit();
	return numRays;
}
#endif


int BvHierarchy::CountViewCells(const ObjectContainer &objects)// const
{
	int nViewCells = 0;

	ViewCell::NewMail();
	ObjectContainer::const_iterator oit, oit_end = objects.end();

	// loop through all object and collect view cell pvs of this node
	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		nViewCells += CountViewCells(*oit);
	}

	return nViewCells;
}


void BvHierarchy::CollectDirtyCandidates(BvhSubdivisionCandidate *sc, 
										 vector<SubdivisionCandidate *> &dirtyList, 
										 const bool onlyUnmailed)
{
	BvhTraversalData &tData = sc->mParentData;
	BvhLeaf *node = tData.mNode;
	
	ViewCellContainer viewCells;
	//ViewCell::NewMail();
	int numRays = CollectViewCells(*tData.mSampledObjects, viewCells, false, false);

	if (0) cout << "collected " << (int)viewCells.size() << " dirty candidates" << endl;
	
	// split candidates handling 
	// these view cells  are thrown into dirty list
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
        VspViewCell *vc = static_cast<VspViewCell *>(*vit);
		VspLeaf *leaf = vc->mLeaves[0];
	
		SubdivisionCandidate *candidate = leaf->GetSubdivisionCandidate();
		
		// is this leaf still a split candidate?
		if (candidate && (!onlyUnmailed || !candidate->Mailed()))
		{
			candidate->Mail();
			candidate->SetDirty(true);
			dirtyList.push_back(candidate);
		}
	}
}


BvhNode *BvHierarchy::GetRoot() const
{
	return mRoot;
}


bool BvHierarchy::IsObjectInLeaf(BvhLeaf *leaf, Intersectable *object) const
{
	ObjectContainer::const_iterator oit =
		lower_bound(leaf->mObjects.begin(), leaf->mObjects.end(), object, ilt);
				
	// objects sorted by id
	if ((oit != leaf->mObjects.end()) && ((*oit)->GetId() == object->GetId()))
	{
		return true;
	}
	else
	{
		return false;
	}
}

#if 0
BvhLeaf *BvHierarchy::GetLeaf(Intersectable *object, BvhNode *node) const
{
	// hack: we use the simpler but faster version
	if (!object) 
		return NULL;

	return object->mBvhLeaf;
	
	///////////////////////////////////////
	// start from root of tree

	if (node == NULL)
		node = mRoot;
	
	vector<BvhLeaf *> leaves;

	stack<BvhNode *> nodeStack;
	nodeStack.push(node);
  
	BvhLeaf *leaf = NULL;
  
	while (!nodeStack.empty())  
	{
		BvhNode *node = nodeStack.top();
		nodeStack.pop();
	
		if (node->IsLeaf()) 
		{
			leaf = static_cast<BvhLeaf *>(node);

			if (IsObjectInLeaf(leaf, object))
			{
				return leaf;
			}
		} 
		else 	
		{	
			// find point
			BvhInterior *interior = static_cast<BvhInterior *>(node);
        
			if (interior->GetBack()->GetBoundingBox().Includes(object->GetBox()))
			{
				nodeStack.push(interior->GetBack());
			}
			
			// search both sides as we are using bounding volumes
			if (interior->GetFront()->GetBoundingBox().Includes(object->GetBox()))
			{
				nodeStack.push(interior->GetFront());
			}
		}
	}
  
	return leaf;
}
#endif

bool BvHierarchy::Export(OUT_STREAM &stream)
{
	ExportNode(mRoot, stream);

	return true;
}


void BvHierarchy::ExportObjects(BvhLeaf *leaf, OUT_STREAM &stream)
{
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		stream << (*oit)->GetId() << " ";
	}
}


void BvHierarchy::ExportNode(BvhNode *node, OUT_STREAM &stream)
{
	if (node->IsLeaf())
	{
		BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
		const AxisAlignedBox3 box = leaf->GetBoundingBox();
		stream << "<Leaf id=\"" << node->GetId() << "\""
			   << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
			   << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\""
			   << " objects=\"";
		
		//-- export objects
		// tmp matt
		if (1) ExportObjects(leaf, stream);
		
		stream << "\" />" << endl;
	}
	else
	{	
		BvhInterior *interior = static_cast<BvhInterior *>(node);
		const AxisAlignedBox3 box = interior->GetBoundingBox();

		stream << "<Interior id=\"" << node->GetId() << "\""
			   << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
			   << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z
			   << "\">" << endl;

		ExportNode(interior->GetBack(), stream);
		ExportNode(interior->GetFront(), stream);

		stream << "</Interior>" << endl;
	}
}


float BvHierarchy::EvalViewCellsVolume(const ObjectContainer &objects)// const
{
	float vol = 0;

	ViewCellContainer viewCells;
	
	// we have to account for all view cells that can 
	// be seen from the objects
	int numRays = CollectViewCells(objects, viewCells, false, false);

	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
		vol += (*vit)->GetVolume();
	}

	return vol;
}


void BvHierarchy::Initialise(const ObjectContainer &objects)
{
	AxisAlignedBox3 box = EvalBoundingBox(objects);

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

	BvhLeaf *bvhleaf = new BvhLeaf(box, NULL, (int)objects.size());
	bvhleaf->mObjects = objects;
	mRoot = bvhleaf;

	// compute bounding box from objects
	mBoundingBox = mRoot->GetBoundingBox();

	// associate root with current objects
	AssociateObjectsWithLeaf(bvhleaf);
}


void BvHierarchy::StoreSampledObjects(ObjectContainer &sampledObjects, const ObjectContainer &objects)
{
	ObjectContainer::const_iterator oit, oit_end = objects.end();

	for (oit = objects.begin(); oit != objects.end(); ++ oit)
	{
		Intersectable *obj = *oit;

		if (!obj->GetOrCreateRays()->empty())
	{
			sampledObjects.push_back(obj);
		}
	}
	}


void BvHierarchy::PrepareConstruction(SplitQueue &tQueue,
									  const VssRayContainer &sampleRays,
									  const ObjectContainer &objects)
{
	///////////////////////////////////////
	//-- we assume that we have objects sorted by their id => 
	//-- we don't have to sort them here and an binary search 
	//-- for identifying if a object is in a leaf.
	
	mBvhStats.Reset();
	mBvhStats.Start();
	mBvhStats.nodes = 1;
		
	// store pointer to this tree
	BvhSubdivisionCandidate::sBvHierarchy = this;
	
	// root and bounding box was already constructed
	BvhLeaf *bvhLeaf = static_cast<BvhLeaf *>(mRoot);
	
	// only rays intersecting objects in node are interesting
	const int nRays = AssociateObjectsWithRays(sampleRays);
	//cout << "using " << nRays << " of " << (int)sampleRays.size() << " rays" << endl;
	
	ObjectContainer *sampledObjects = new ObjectContainer();
	StoreSampledObjects(*sampledObjects, objects);

#if STORE_VIEWCELLS_WITH_BVH
	AssociateViewCellsWithObjects(*sampledObjects);
#endif

	// probability that volume is "seen" from the view cells
	const float prop = EvalViewCellsVolume(*sampledObjects) / GetViewSpaceVolume();

	// create bvh traversal data
	BvhTraversalData oData(bvhLeaf, 0, prop, nRays);
		
	// create sorted object lists for the first data
	if (mUseGlobalSorting)
	{
		AssignInitialSortedObjectList(oData, objects);
	}
	
	oData.mSampledObjects = sampledObjects;
	
	///////////////////
	//-- add first candidate for object space partition	

	mTotalCost = EvalRenderCost(objects);
	mPvsEntries = CountViewCells(*sampledObjects);

	oData.mCorrectedPvs = oData.mPvs = (float)mPvsEntries;
	oData.mCorrectedVolume = oData.mVolume = prop;
	
	BvhSubdivisionCandidate *oSubdivisionCandidate = 
		new BvhSubdivisionCandidate(oData);

	bvhLeaf->SetSubdivisionCandidate(oSubdivisionCandidate);

#if STORE_VIEWCELLS_WITH_BVH
	ReleaseViewCells(*sampledObjects);
#endif

	if (mApplyInitialPartition)
	{
		vector<SubdivisionCandidate *> candidateContainer;

		mIsInitialSubdivision = true;
		
		// evaluate priority
		EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
		PrintSubdivisionStats(*oSubdivisionCandidate);

		ApplyInitialSubdivision(oSubdivisionCandidate, candidateContainer);		

		mIsInitialSubdivision = false;

		vector<SubdivisionCandidate *>::const_iterator cit, cit_end = candidateContainer.end();

		for (cit = candidateContainer.begin(); cit != cit_end; ++ cit)
		{
			BvhSubdivisionCandidate *sCandidate = static_cast<BvhSubdivisionCandidate *>(*cit);
			
			// reevaluate priority
			EvalSubdivisionCandidate(*sCandidate, true, true);
			tQueue.Push(sCandidate);
		}

		cout << "size of initial bv subdivision: " << GetStatistics().Leaves() << endl;
	}
	else
	{	
		// evaluate priority
		EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
		PrintSubdivisionStats(*oSubdivisionCandidate);

		tQueue.Push(oSubdivisionCandidate);
		cout << "size of initial bv subdivision: " << GetStatistics().Leaves() << endl;
	}
}


void BvHierarchy::AssignInitialSortedObjectList(BvhTraversalData &tData,
												const ObjectContainer &objects)
{
	const bool doSort = true;

	// we sort the objects as a preprocess so they don't have
	// to be sorted for each split
	for (int i = 0; i < 3; ++ i)
	{
		SortableEntryContainer *sortedObjects = new SortableEntryContainer();

		CreateLocalSubdivisionCandidates(objects, 
									     &sortedObjects, 
										 doSort, 
										 i);
		
		// copy list into traversal data list
		tData.mSortedObjects[i] = new ObjectContainer();
		tData.mSortedObjects[i]->reserve((int)objects.size());

		SortableEntryContainer::const_iterator oit, oit_end = sortedObjects->end();

		for (oit = sortedObjects->begin(); oit != oit_end; ++ oit)
		{
			tData.mSortedObjects[i]->push_back((*oit).mObject);
		}

		delete sortedObjects;
	}

	// next sorted list: by size (for initial criteria)
	tData.mSortedObjects[3] = new ObjectContainer();
	tData.mSortedObjects[3]->reserve((int)objects.size());

	*(tData.mSortedObjects[3]) = objects;
	
	stable_sort(tData.mSortedObjects[3]->begin(), tData.mSortedObjects[3]->end(), smallerSize);
}


void BvHierarchy::AssignSortedObjects(const BvhSubdivisionCandidate &sc,
									  BvhTraversalData &frontData,
									  BvhTraversalData &backData)
{
	Intersectable::NewMail();

	// we sorted the objects as a preprocess so they don't have
	// to be sorted for each split
	ObjectContainer::const_iterator fit, fit_end = sc.mFrontObjects.end();

	for (fit = sc.mFrontObjects.begin(); fit != fit_end; ++ fit)
	{
		(*fit)->Mail();
	}

	for (int i = 0; i < 4; ++ i)
	{
		frontData.mSortedObjects[i] = new ObjectContainer();
		backData.mSortedObjects[i] = new ObjectContainer();

		frontData.mSortedObjects[i]->reserve((int)sc.mFrontObjects.size());
		backData.mSortedObjects[i]->reserve((int)sc.mBackObjects.size());

		ObjectContainer::const_iterator oit, oit_end = sc.mParentData.mSortedObjects[i]->end();

		// all the front objects are mailed => assign the sorted object lists
		for (oit = sc.mParentData.mSortedObjects[i]->begin(); oit != oit_end; ++ oit)
		{
			if ((*oit)->Mailed())
			{
				frontData.mSortedObjects[i]->push_back(*oit);
			}
			else
			{
				backData.mSortedObjects[i]->push_back(*oit);
			}
		}
	}
}


void BvHierarchy::Reset(SplitQueue &tQueue,
						const VssRayContainer &sampleRays, 
						const ObjectContainer &objects)
{

	// reset stats
	mBvhStats.Reset();
	mBvhStats.Start();
	mBvhStats.nodes = 1;

	// reset root
	DEL_PTR(mRoot);
	
	BvhLeaf *bvhleaf = new BvhLeaf(mBoundingBox, NULL, (int)objects.size());
	bvhleaf->mObjects = objects;
	mRoot = bvhleaf;
	
	ObjectContainer *sampledObjects = new ObjectContainer();
	StoreSampledObjects(*sampledObjects, objects);

#if STORE_VIEWCELLS_WITH_BVH
	AssociateViewCellsWithObjects(*sampledObjects);
#endif

	//mTermMinProbability *= mVspTree->GetBoundingBox().GetVolume();
	// probability that volume is "seen" from the view cells
	const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
	const float prop = EvalViewCellsVolume(*sampledObjects);

	const int nRays = CountRays(*sampledObjects);
	BvhLeaf *bvhLeaf = static_cast<BvhLeaf *>(mRoot);

	// create bvh traversal data
	BvhTraversalData oData(bvhLeaf, 0, prop, nRays);

	oData.mSampledObjects = sampledObjects;

	if (mUseGlobalSorting)
		AssignInitialSortedObjectList(oData, objects);
	
#if STORE_VIEWCELLS_WITH_BVH
	ReleaseViewCells(*sampledObjects);
#endif
	///////////////////
	//-- add first candidate for object space partition	

	BvhSubdivisionCandidate *oSubdivisionCandidate = 
		new BvhSubdivisionCandidate(oData);

	EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
	bvhLeaf->SetSubdivisionCandidate(oSubdivisionCandidate);

	mTotalCost = (float)objects.size() * prop;

	PrintSubdivisionStats(*oSubdivisionCandidate);

	tQueue.Push(oSubdivisionCandidate);
}


void BvhStatistics::Print(ostream &app) const
{
	app << "=========== BvHierarchy statistics ===============\n";

	app << setprecision(4);

	app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";

	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";

	app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";

	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";

	app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits << endl;

	app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
		<< maxCostNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
		<< minProbabilityNodes * 100 / (double)Leaves() << endl;


	//////////////////////////////////////////////////
	
	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
		<<	maxDepthNodes * 100 / (double)Leaves() << endl;
	
	app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;

	app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;

	app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;

	
	////////////////////////////////////////////////////////
	
	app << "#N_PMINOBJECTSLEAVES  ( Percentage of leaves with mininum objects )\n"
		<< minObjectsNodes * 100 / (double)Leaves() << endl;

	app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" << maxObjectRefs << "\n";

	app << "#N_MINOBJECTREFS  ( Min number of object refs / leaf )\n" << minObjectRefs << "\n";

	app << "#N_EMPTYLEAFS ( Empty leafs )\n" << emptyNodes << "\n";
	
	app << "#N_PAVGOBJECTSLEAVES  ( average object refs / leaf)\n" << AvgObjectRefs() << endl;


	////////////////////////////////////////////////////////
	
	app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with mininum rays )\n"
		<< minRaysNodes * 100 / (double)Leaves() << endl;

	app << "#N_MAXRAYREFS  ( Max number of ray refs / leaf )\n" << maxRayRefs << "\n";

	app << "#N_MINRAYREFS  ( Min number of ray refs / leaf )\n" << minRayRefs << "\n";
	
	app << "#N_PAVGRAYLEAVES  ( average ray refs / leaf )\n" << AvgRayRefs() << endl;
	
	app << "#N_PAVGRAYCONTRIBLEAVES  ( Average ray contribution)\n" <<
		rayRefs / (double)objectRefs << endl;

	app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
		maxRayContriNodes * 100 / (double)Leaves() << endl;

	app << "#N_PGLOBALCOSTMISSES ( Global cost misses )\n" << mGlobalCostMisses << endl;

	app << "========== END OF BvHierarchy statistics ==========\n";
}


// TODO: return memory usage in MB
float BvHierarchy::GetMemUsage() const
{
	return (float)(sizeof(BvHierarchy)
				   + mBvhStats.Leaves() * sizeof(BvhLeaf) 
				   + mBvhStats.Interior() * sizeof(BvhInterior)
				   ) / float(1024 * 1024);
}


void BvHierarchy::SetActive(BvhNode *node) const
{
	vector<BvhLeaf *> leaves;

	// sets the pointers to the currently active view cells
	CollectLeaves(node, leaves);
	vector<BvhLeaf *>::const_iterator lit, lit_end = leaves.end();

	for (lit = leaves.begin(); lit != lit_end; ++ lit)
	{
		(*lit)->SetActiveNode(node);
	}
}


BvhNode *BvHierarchy::SubdivideAndCopy(SplitQueue &tQueue, 
									   SubdivisionCandidate *splitCandidate)
{
	BvhSubdivisionCandidate *sc = 
		static_cast<BvhSubdivisionCandidate *>(splitCandidate);
	BvhTraversalData &tData = sc->mParentData;

	BvhNode *currentNode = tData.mNode;
	BvhNode *oldNode = (BvhNode *)splitCandidate->mEvaluationHack;

	if (!oldNode->IsLeaf())
	{	
		//////////////
		//-- continue subdivision

		BvhTraversalData tFrontData;
		BvhTraversalData tBackData;
			
		BvhInterior *oldInterior = static_cast<BvhInterior *>(oldNode);
		
		sc->mFrontObjects.clear();
		sc->mBackObjects.clear();

		oldInterior->GetFront()->CollectObjects(sc->mFrontObjects);
		oldInterior->GetBack()->CollectObjects(sc->mBackObjects);
		
		// evaluate the changes in render cost and pvs entries
		EvalSubdivisionCandidate(*sc, false, true);

		// create new interior node and two leaf node
		currentNode = SubdivideNode(*sc, tFrontData, tBackData);
	
		//oldNode->mRenderCostDecr += sc->GetRenderCostDecrease();
		//oldNode->mPvsEntriesIncr += sc->GetPvsEntriesIncr();
		
		//oldNode->mRenderCostDecr = sc->GetRenderCostDecrease();
		//oldNode->mPvsEntriesIncr = sc->GetPvsEntriesIncr();
		
		///////////////////////////
		//-- push the new split candidates on the queue
		
		BvhSubdivisionCandidate *frontCandidate = new BvhSubdivisionCandidate(tFrontData);
		BvhSubdivisionCandidate *backCandidate = new BvhSubdivisionCandidate(tBackData);

		frontCandidate->SetPriority((float)-oldInterior->GetFront()->GetTimeStamp());
		backCandidate->SetPriority((float)-oldInterior->GetBack()->GetTimeStamp());

		frontCandidate->mEvaluationHack = oldInterior->GetFront();
		backCandidate->mEvaluationHack = oldInterior->GetBack();

		// cross reference
		tFrontData.mNode->SetSubdivisionCandidate(frontCandidate); 
		tBackData.mNode->SetSubdivisionCandidate(backCandidate);

		//cout << "f: " << frontCandidate->GetPriority() << " b: " << backCandidate->GetPriority() << endl;
		tQueue.Push(frontCandidate);
		tQueue.Push(backCandidate);
	}

	/////////////////////////////////
	//-- node is a leaf => terminate traversal

	if (currentNode->IsLeaf())
	{
		// this leaf is no candidate for splitting anymore
		// => detach subdivision candidate
		tData.mNode->SetSubdivisionCandidate(NULL); 
		// detach node so we don't delete it with the traversal data
		tData.mNode = NULL;
	}
	
	return currentNode;
}


void BvHierarchy::CollectObjects(const AxisAlignedBox3 &box,
								 ObjectContainer &objects)
{
  stack<BvhNode *> nodeStack;
  
  nodeStack.push(mRoot);

  while (!nodeStack.empty()) {
	BvhNode *node = nodeStack.top();
	
	nodeStack.pop();
	if (node->IsLeaf()) {
	  BvhLeaf *leaf = (BvhLeaf *)node;
	  if (Overlap(box, leaf->GetBoundingBox())) {
		Intersectable *object = leaf;
		if (!object->Mailed()) {
		  object->Mail();
		  objects.push_back(object);
		}
	  }
	} 
	else 
	  {
		BvhInterior *interior = (BvhInterior *)node;
		if (Overlap(box, interior->GetBoundingBox())) {
		  bool pushed = false;
		  if (!interior->GetFront()->Mailed()) {
			nodeStack.push(interior->GetFront());
			pushed = true;
		  }
		  if (!interior->GetBack()->Mailed()) {
			nodeStack.push(interior->GetBack());
			pushed = true;
		  }
		  // avoid traversal of this node in the next query
		  if (!pushed)
			interior->Mail();
		}
	  }
  }
}


void BvHierarchy::CreateUniqueObjectIds()
{
	stack<BvhNode *> nodeStack;
	nodeStack.push(mRoot);

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

		node->SetId(currentId ++);

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

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


void BvHierarchy::ApplyInitialSubdivision(SubdivisionCandidate *firstCandidate,
										  vector<SubdivisionCandidate *> &candidateContainer)
{
	SplitQueue tempQueue;
	tempQueue.Push(firstCandidate);

	while (!tempQueue.Empty())
	{
		SubdivisionCandidate *candidate = tempQueue.Top();
		tempQueue.Pop();

		BvhSubdivisionCandidate *bsc = 
			static_cast<BvhSubdivisionCandidate *>(candidate);

		if (!InitialTerminationCriteriaMet(bsc->mParentData))
		{
			const bool globalCriteriaMet = GlobalTerminationCriteriaMet(bsc->mParentData);
		
			SubdivisionCandidateContainer dirtyList;
			BvhNode *node = Subdivide(tempQueue, bsc, globalCriteriaMet, dirtyList);

			// not needed anymore
			delete bsc;
		}
		else 
		{
			// initial preprocessing  finished for this candidate
			// add to candidate container
			candidateContainer.push_back(bsc);
		}
	}
}


void BvHierarchy::ApplyInitialSplit(const BvhTraversalData &tData,
									ObjectContainer &frontObjects,
									ObjectContainer &backObjects)
{
	ObjectContainer *objects = tData.mSortedObjects[3];

	ObjectContainer::const_iterator oit, oit_end = objects->end();
   
	float maxAreaDiff = -1.0f;

	ObjectContainer::const_iterator backObjectsStart = objects->begin();

	for (oit = objects->begin(); oit != (objects->end() - 1); ++ oit)
	{
		Intersectable *objS = *oit;
		Intersectable *objL = *(oit + 1);
		
		const float areaDiff = 
				objL->GetBox().SurfaceArea() - objS->GetBox().SurfaceArea();

		if (areaDiff > maxAreaDiff)
		{
			maxAreaDiff = areaDiff;
			backObjectsStart = oit + 1;
		}
	}

	// belongs to back bv
	for (oit = objects->begin(); oit != backObjectsStart; ++ oit)
	{
		frontObjects.push_back(*oit);
	}

	// belongs to front bv
	for (oit = backObjectsStart; oit != oit_end; ++ oit)
	{
		backObjects.push_back(*oit);
	}
	
	cout << "front: " << (int)frontObjects.size() << " back: " << (int)backObjects.size() << " " 
		 << backObjects.front()->GetBox().SurfaceArea() - frontObjects.back()->GetBox().SurfaceArea() << endl;
}


inline static float AreaRatio(Intersectable *smallObj, Intersectable *largeObj)
{
	const float areaSmall = smallObj->GetBox().SurfaceArea();
	const float areaLarge = largeObj->GetBox().SurfaceArea();

	return areaSmall / (areaLarge - areaSmall + Limits::Small);
}


bool BvHierarchy::InitialTerminationCriteriaMet(const BvhTraversalData &tData) const
{
	const bool terminationCriteriaMet =
			(0
		    || ((int)tData.mNode->mObjects.size() < mInitialMinObjects)
			|| (tData.mNode->mObjects.back()->GetBox().SurfaceArea() < mInitialMinArea)
			|| (AreaRatio(tData.mNode->mObjects.front(), tData.mNode->mObjects.back()) > mInitialMaxAreaRatio)
			);

	cout << "criteria met: "<< terminationCriteriaMet << "\n" 
		 << "size: " << (int)tData.mNode->mObjects.size() << " max: " << mInitialMinObjects << endl
		 << "ratio: " << AreaRatio(tData.mNode->mObjects.front(), tData.mNode->mObjects.back()) << " max: " << mInitialMaxAreaRatio << endl
		 << "area: " << tData.mNode->mObjects.back()->GetBox().SurfaceArea() << " max: " << mInitialMinArea << endl << endl;

	return terminationCriteriaMet;
}


// HACK
float BvHierarchy::GetRenderCostIncrementially(BvhNode *node) const
{
	if (node->mRenderCost < 0)
	{
		//cout <<"p";
		if (node->IsLeaf())
		{
			BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
			node->mRenderCost = EvalAbsCost(leaf->mObjects);
		}
		else
		{
			BvhInterior *interior = static_cast<BvhInterior *>(node);
		
			node->mRenderCost = GetRenderCostIncrementially(interior->GetFront()) + 
								GetRenderCostIncrementially(interior->GetBack());
		}
	}

	return node->mRenderCost;
}


void BvHierarchy::Compress()
{
}


void BvHierarchy::SetUniqueNodeIds()
{
	// export bounding boxes
	vector<BvhNode *> nodes;

	// hack: should also expect interior nodes
	CollectNodes(mRoot, nodes);

	vector<BvhNode *>::const_iterator oit, oit_end = nodes.end();

	int id = 0;

	for (oit = nodes.begin(); oit != oit_end; ++ oit, ++ id)
	{
		(*oit)->SetId(id);
	}
}


}