#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"


namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0

static float debugVol = 0;
int BvhNode::sMailId = 0;
BvHierarchy *BvHierarchy::BvhSubdivisionCandidate::sBvHierarchy = NULL;


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


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

BvhNode::BvhNode(): mParent(NULL)
{
}

BvhNode::BvhNode(const AxisAlignedBox3 &bbox):
mBoundingBox(bbox)
{
}


BvhNode::BvhNode(const AxisAlignedBox3 &bbox, BvhInterior *parent):
mBoundingBox(bbox), mParent(parent)
{
}


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


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


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



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


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox):
BvhNode(bbox)
{
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent):
BvhNode(bbox, parent)
{
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent, const int numObjects):
BvhNode(bbox, parent)
{
	mObjects.reserve(numObjects);
}


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


BvhLeaf::~BvhLeaf() 
{ 
}


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


BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox):
BvhNode(bbox)
{
}


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


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;
}



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


BvHierarchy::BvHierarchy():
mRoot(NULL),
mTimeStamp(1)
{
	ReadEnvironment();
	mSubdivisionCandidates = new vector<SortableEntry>;
}


BvHierarchy::~BvHierarchy()
{
	// delete kd intersectables
	BvhIntersectableMap::iterator it, it_end = mBvhIntersectables.end();

	for (it = mBvhIntersectables.begin(); it != mBvhIntersectables.end(); ++ it)
	{
		DEL_PTR((*it).second);
	}

	DEL_PTR(mSubdivisionCandidates);

	mSubdivisionStats.close();
}


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

	//-- termination criteria for autopartition
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxLeaves", mTermMaxLeaves);
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.minObjects", mTermMinObjects);
	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minProbability", mTermMinVolume);
	
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.missTolerance", mTermMissTolerance);
	
	//-- max cost ratio for early tree termination
	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.maxCostRatio", mTermMaxCostRatio);

	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minGlobalCostRatio",
		mTermMinGlobalCostRatio);
	

	//-- factors for bsp tree split plane heuristics
	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.ct_div_ci", mCtDivCi);

	//-- partition criteria
	Environment::GetSingleton()->GetFloatValue("OspTree.Construction.epsilon", mEpsilon);

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


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

	Environment::GetSingleton()->GetFloatValue("OspTree.Construction.splitBorder", mSplitBorder);
	Environment::GetSingleton()->GetFloatValue("OspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
	

	//-- debug output

	Debug << "******* Bvh hierarchy options ******** " << endl;

    Debug << "max depth: " << mTermMaxDepth << endl;
	Debug << "min probabiliy: " << mTermMinVolume<< 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 << "subdivision stats log: " << subdivisionStatsLog << endl;
	
	Debug << "split borders: " << mSplitBorder << endl;
	Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
	Debug << endl;
}


void AssociateObjectsWithLeaf(BvhLeaf *leaf)
{
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		(*oit)->mBvhLeaf = leaf;
	}
}


BvhInterior *BvHierarchy::SubdivideNode(const ObjectContainer &frontObjects,
										const ObjectContainer &backObjects,
										const BvhTraversalData &tData,
										BvhTraversalData &frontData,
										BvhTraversalData &backData)
{
	BvhLeaf *leaf = tData.mNode;
	
	// two new leaves
    mBvhStats.nodes += 2;

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


	//-- the front and back traversal data is filled with the new values

	frontData.mDepth = backData.mDepth = tData.mDepth + 1;
		
	// frontData.mRays = new RayInfoContainer();
	// backData.mRays = new RayInfoContainer();

	frontData.mBoundingBox = ComputeBoundingBox(frontObjects);
	backData.mBoundingBox = ComputeBoundingBox(backObjects);
	//RemoveParentViewCellsPvs(leaf, *tData.mRays);


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

	BvhLeaf *back = new BvhLeaf(backData.mBoundingBox, node, (int)backObjects.size());
	BvhLeaf *front = new BvhLeaf(frontData.mBoundingBox, node, (int)frontObjects.size());

	BvhInterior *parent = leaf->GetParent();


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

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

	++ mBvhStats.splits;

	//UpdateViewCellsPvs(front, *frontData.mRays);
	//UpdateViewCellsPvs(back, *backData.mRays);

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

	//ProcessMultipleRefs(front);
    //ProcessMultipleRefs(back);

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

	back->mObjects = backObjects;
	front->mObjects = frontObjects;

	AssociateObjectsWithLeaf(back);
	AssociateObjectsWithLeaf(front);
   

	// compute probability, i.e., volume of seen view cells
	frontData.mVolume = frontData.mBoundingBox.GetVolume();
	backData.mVolume =  backData.mBoundingBox.GetVolume();

	return node;
}


BvhNode *BvHierarchy::Subdivide(SplitQueue &tQueue,
								SubdivisionCandidate *splitCandidate,
								const bool globalCriteriaMet)
{
	BvhSubdivisionCandidate *sc = dynamic_cast<BvhSubdivisionCandidate *>(splitCandidate);
	BvhTraversalData &tData = sc->mParentData;

	BvhNode *newNode = tData.mNode;

	if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
	{	
		BvhTraversalData tFrontData;
		BvhTraversalData tBackData;

		//-- continue subdivision
		
		// create new interior node and two leaf node
		newNode = SubdivideNode(sc->mFrontObjects,
								sc->mBackObjects,
								tData, 
								tFrontData, 
								tBackData);
	
		const int maxCostMisses = sc->mMaxCostMisses;

        // how often was max cost ratio missed in this branch?
		tFrontData.mMaxCostMisses = maxCostMisses;
		tBackData.mMaxCostMisses = maxCostMisses;
		
		mTotalCost -= sc->GetRenderCostDecrease();

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

		//-- push the new split candidates on the queue
		
		BvhSubdivisionCandidate *frontCandidate = new BvhSubdivisionCandidate(tFrontData);
		BvhSubdivisionCandidate *backCandidate = new BvhSubdivisionCandidate(tBackData);

		EvalSubdivisionCandidate(*frontCandidate);
		EvalSubdivisionCandidate(*backCandidate);

		tQueue.Push(frontCandidate);
		tQueue.Push(backCandidate);

		// delete old leaf node
		DEL_PTR(tData.mNode);
	}


	//-- terminate traversal
	if (newNode->IsLeaf())
	{
		EvaluateLeafStats(tData);
        
		const bool mStoreRays= true;

		//-- store additional info
		if (mStoreRays)
		{
			BvhLeaf *leaf = dynamic_cast<BvhLeaf *>(newNode);
			CollectRays(leaf->mObjects, leaf->mVssRays);
		}
	}
	
	//-- cleanup
	tData.Clear();
	
	return newNode;
}


void BvHierarchy::EvalSubdivisionCandidate(BvhSubdivisionCandidate &splitCandidate)
{
	ObjectContainer frontObjects, backObjects;

	// compute best object partition
	const float ratio =
		SelectObjectPartition(splitCandidate.mParentData, frontObjects, backObjects);

	const bool success = ratio < mTermMaxCostRatio;

	float oldRenderCost;

	// compute global decrease in render cost
	const float renderCostDecr = EvalRenderCostDecrease(frontObjects, 
														backObjects,
														splitCandidate.mParentData,
														oldRenderCost);

	splitCandidate.SetRenderCostDecrease(renderCostDecr);

#if 0
	const float priority = (float)-data.mDepth;
#else	
	// take render cost of node into account 
	// otherwise danger of being stuck in a local minimum!!
	const float factor = mRenderCostDecreaseWeight;
	const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;

#endif

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

	splitCandidate.mMaxCostMisses = 
		success ? splitCandidate.mParentData.mMaxCostMisses : 
	splitCandidate.mParentData.mMaxCostMisses + 1;
}


inline bool BvHierarchy::LocalTerminationCriteriaMet(const BvhTraversalData &data) const
{
	// matt: TODO
	return ( 
		//(data.mNode->mObjects.size() < mTermMinObjects) ||
		//(data.mProbability <= mTermMinProbability) ||
		(data.mDepth >= mTermMaxDepth)
		 );
}


inline bool BvHierarchy::GlobalTerminationCriteriaMet(const BvhTraversalData &data) const
{
	// matt: TODO
	return (
		(mBvhStats.Leaves() >= mTermMaxLeaves)
		//mOutOfMemory || 
		//(mGlobalCostMisses >= mTermGlobalCostMissTolerance)
		);
}


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

	if (data.mDepth >= mTermMaxDepth)
	{
        ++ mBvhStats.maxDepthNodes;
		//Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
	}

	if (data.mVolume <= mTermMinVolume)
		++ mBvhStats.minProbabilityNodes;
	
	// accumulate depth to compute average depth
	mBvhStats.accumDepth += data.mDepth;
 
	//	if ((int)(leaf->mObjects.size()) < mTermMinCost)
	//     ++ mOspStats.minCostNodes;
 
	if ((int)(leaf->mObjects.size()) > mBvhStats.maxObjectRefs)
		mBvhStats.maxObjectRefs = (int)leaf->mObjects.size();
}


float BvHierarchy::EvalLocalCostHeuristics(const BvhTraversalData &tData,
										   const int axis,
										   ObjectContainer &objectsFront,
										   ObjectContainer &objectsBack)
{
	// the relative cost ratio
	float ratio = 99999999.0f;

	// 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
	
	const float minBox = tData.mBoundingBox.Min(axis);
	const float maxBox = tData.mBoundingBox.Max(axis);

	const float sizeBox = maxBox - minBox;

	float minBand = minBox + mSplitBorder * (maxBox - minBox);
	float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);

	//-- sort so we can use a sweep
	SortSubdivisionCandidates(tData, axis, minBand, maxBand);

	float totalVol = PrepareHeuristics(tData);
	float voll = 0;
	float volr = totalVol;

	/// this is kind of a reverse pvs
	const int totalObjects = (int)tData.mNode->mObjects.size();
	
	int objectsl = 0;
	int objectsr = totalObjects;

	float sum = (float)totalVol * objectsr;


	bool splitPlaneFound = false;

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

	float volBack = voll;
	float volFront = volr;

	int nObjectsBack = objectsl;
	int nObjectsFront = objectsr;
	
	float minSum = 1e20f;

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

	debugVol = 0;
	const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();

	/////////////////////////////
	// the sweep heuristics

	// mail the objects on the left side
	Intersectable::NewMail();
	
	//-- traverse through events and find best split plane
	vector<SortableEntry>::const_iterator ci, ci_end = mSubdivisionCandidates->end();

	for (ci = mSubdivisionCandidates->begin(); ci != ci_end; ++ ci)
	{
		Intersectable *object = (*ci).mObject;

		EvalHeuristicsContribution(tData.mNode, *ci, voll, volr, objectsl);
		objectsr = totalObjects - objectsl;

		// Note: sufficient to compare size of bounding boxes of front and back side?

		if (((*ci).mPos >= minBand) && ((*ci).mPos <= maxBand))
		{
			// voll = view cells that see only left node (i.e., left pvs)
			// volr = view cells that see only right node (i.e., right pvs)
			// rest = view cells that see both nodes (i.e., total pvs)
            sum = voll * objectsl + volr * objectsr + (totalVol - voll - volr) * totalObjects;

			float currentPos;
			
			// HACK: current positition is BETWEEN visibility events
			if (1 && ((ci + 1) != ci_end))
				currentPos = ((*ci).mPos + (*(ci + 1)).mPos) * 0.5f;
			else
				currentPos = (*ci).mPos;			

			if ((totalVol - voll - volr - debugVol) / viewSpaceVol > 0.0001)
				Debug << "front and back volume: " << (totalVol - voll - volr) / viewSpaceVol << " error: " << (totalVol - voll - volr - debugVol) / viewSpaceVol << endl;
#if 0				
			Debug << "pos: " << currentPos 
				 << "\t (pvsl: " << pvsl << ", pvsr: " << pvsr << ")"
				 << "\t (voll: " << voll << ", volr: " << volr << ")"
				 << "\t (vcl: " << vcl << ", vcr: " << vcr << ", nvc: " << numViewCells << ")"
				 << "\t sum: " << sum << endl;
				
#endif
			if (sum < minSum)
			{
				splitPlaneFound = true;

				minSum = sum;

				nObjectsBack = objectsl;
				nObjectsFront = objectsr;
				
				volBack = voll;
				volFront = volr;
			}
		}
	}
	
	//-- compute cost
	const float frontAndBackVol = totalVol - volFront - volBack;

	const float oldRenderCost = (float)totalObjects * totalVol + Limits::Small;
	const float newRenderCost = (float)nObjectsFront * volFront + 
								(float)nObjectsBack * volBack +
								(float)totalObjects * frontAndBackVol;

	if (splitPlaneFound)
	{
		ratio = newRenderCost / oldRenderCost;
	}

	//Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" 
	//	  << position << " t=" << (position - minBox) / (maxBox - minBox)
	//      << "\t pb=(" << volBack << ")\t pf=(" << volFront << ")" << endl;

	Debug << "\n§§§§ eval local cost §§§§" << endl
		  << "back pvs: " << nObjectsBack << " front pvs: " << nObjectsFront << " total pvs: " << totalObjects << 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;

	if (oldRenderCost < newRenderCost)
		Debug << "\nwarning!!:\n" << "old rc: " << oldRenderCost * viewSpaceVol << " new rc: " << newRenderCost * viewSpaceVol << endl;

	// assign objects
	ObjectContainer::const_iterator oit, oit_end = tData.mNode->mObjects.end();

	for (oit = tData.mNode->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

		if (obj->Mailed()) // belongs to back objects
			objectsBack.push_back(obj);
		else
			objectsFront.push_back(obj);
	}

	return ratio;
}


void BvHierarchy::SortSubdivisionCandidates(const BvhTraversalData &tData,
									  const int axis, 
									  float minBand, 
									  float maxBand)

{
	mSubdivisionCandidates->clear();

	//RayInfoContainer *rays = tData.mRays;
	BvhLeaf *leaf = tData.mNode;

	// compute requested size
	int requestedSize = 0;
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
	for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
		requestedSize += 2 + (int)(*oit)->mVssRays.size() * 2;
	
	// creates a sorted split candidates array
	if (mSubdivisionCandidates->capacity() > 500000 &&
		requestedSize < (int)(mSubdivisionCandidates->capacity() / 10) )
	{
    	delete mSubdivisionCandidates;
   		mSubdivisionCandidates = new vector<SortableEntry>;
	}

	mSubdivisionCandidates->reserve(requestedSize);

	float pos;

	//-- insert object queries
	//-- These queries can induce a change in pvs size.
	//-- Note that they cannot induce a change in view cell volume, as
	//-- there is no ray associated with these events!!

	for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		
		Intersectable *object = *oit;
		const AxisAlignedBox3 box = object->GetBox();

		mSubdivisionCandidates->push_back(
			SortableEntry(SortableEntry::OBJECT,
						  box.Min(axis),
						  object,
						  NULL)
						  );


		//-- insert ray queries
		//-- we are intersested only in rays which intersect an object that
		//-- is part of the kd node because they can induce a change in view cell
		//-- volume on the left and the right part. 
		//-- Note that they cannot induce a change in pvs size!!
		
		VssRayContainer::const_iterator rit, rit_end = obj->mVssRays.end();

		for (rit = obj->mVssRays.begin(); rit < rit_end; ++ rit)
		{
			VssRay *ray = (*rit);
			const bool positive = ray->HasPosDir(axis);
	
			pos = ray->mTermination[axis];

			mSubdivisionCandidates->push_back(
				SortableEntry(SortableEntry::VIEWCELLS, 
				pos, 
				ray->mTerminationObject, 
				ray)
				);
		}
	}

	stable_sort(mSubdivisionCandidates->begin(), mSubdivisionCandidates->end());
}


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


void BvHierarchy::EvalRayContribution(const VssRay &ray, 
									  float &volLeft, 
									  float &volRight)
{
	ViewCellContainer viewCells;
	mVspTree->GetViewCells(ray, viewCells);

	/// 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;

		const float vol = viewCell->GetVolume();

		if (!viewCell->Mailed())
		{
			viewCell->Mail();

			// we now see view cell from both nodes 
			// => remove ref to right child node
			volRight -= vol;
			debugVol += vol;
		}

		// last reference into the right node
		if (-- viewCell->mCounter == 0)
		{
			// view cell was previously seen from both nodes  =>
			// contributes only to left node now
			volLeft += vol;
			debugVol -= vol;
		}
	}
}


float BvHierarchy::PrepareHeuristics(const VssRay &ray)
{
	float vol = 0;

	ViewCellContainer viewCells;
	mVspTree->GetViewCells(ray, viewCells);
    	
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

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

		if (!vc->Mailed())
		{
			//Debug << "single vol: "<< vc->GetVolume() << endl;
			vc->Mail();
			vc->mCounter = 0;
			vol += vc->GetVolume();
		}
		
		// view cell volume already added => just increase counter
		++ vc->mCounter;
	}

	return vol;
}


float BvHierarchy::PrepareHeuristics(const BvhTraversalData &tData)
{	
	float vol = 0;
	debugVol = 0;
	
	ViewCell::NewMail();

	BvhLeaf *leaf = tData.mNode;

	VssRayContainer rays;
	CollectRays(leaf->mObjects, rays);

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

	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);

		const float volContri = PrepareHeuristics(*ray);

		// if hitpoint with one of the objects is inside this node, we
		// evaluate the volume of the view cells seen by this ray
		if (ray->mTerminationObject)
		{
            vol += volContri;
		}

		// count double if both hit points are within the kd node
		if (0 && ray->mOriginObject)
		{
			vol += volContri;
		}
	}

	return vol;
}
///////////////////////////////////////////////////////////


void BvHierarchy::EvalHeuristicsContribution(BvhLeaf *leaf,
											 const SortableEntry &ci,
											 float &volLeft, 
											 float &volRight, 
											 int &objectsLeft)
{
	Intersectable *obj = ci.mObject;
	VssRay *ray = ci.mRay;

	// switch between different types of events
	switch (ci.mType) 
	{
		case SortableEntry::OBJECT:
			cout << "o";
			ci.mObject->Mail();
			++ objectsLeft;
			break;

		// compute volume contribution from view cells
		case SortableEntry::VIEWCELLS:
			cout << "v";
			// process ray if the hit point with termination / origin object 
			// is inside this kd leaf
			EvalRayContribution(*ray, volLeft, volRight);
			break;

		default:
			cout << "should not come here" << endl;
			break;
	}
	//cout << "vol left: " << volLeft << " vol right " << volRight << endl;
}


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


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


float BvHierarchy::SelectObjectPartition(const BvhTraversalData &tData,
										 ObjectContainer &frontObjects,
										 ObjectContainer &backObjects)
{
	ObjectContainer nFrontObjects[3];
	ObjectContainer nBackObjects[3];

	float nCostRatio[3];

	// create bounding box of node geometry
	AxisAlignedBox3 box = tData.mBoundingBox;
		
	int sAxis = 0;
	int bestAxis = -1;

	if (mOnlyDrivingAxis)
	{
		sAxis = box.Size().DrivingAxis();
	}
	
	// -- evaluate split cost for all three axis
	
	for (int axis = 0; axis < 3; ++ axis)
	{
		if (!mOnlyDrivingAxis || (axis == sAxis))
		{
			if (1 || mUseCostHeuristics)
			{
				//-- partition objects using heuristics
				
				nCostRatio[axis] =
					EvalLocalCostHeuristics(
											tData,
											axis,
											nFrontObjects[axis],
											nBackObjects[axis]);			
			}
			if (bestAxis == -1)
			{
				bestAxis = axis;
			}
			else if (nCostRatio[axis] < nCostRatio[bestAxis])
			{
				bestAxis = axis;
			} 
		}
	}

	//-- assign values
	
	frontObjects = nFrontObjects[bestAxis];
	backObjects = nFrontObjects[bestAxis];

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


void BvHierarchy::AssociateObjectsWithRays(const VssRayContainer &rays)
{

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

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

		if (ray->mTerminationObject)
		{
			ray->mTerminationObject->mVssRays.push_back(ray);
		}

		if (0 && ray->mOriginObject)
		{
			ray->mOriginObject->mVssRays.push_back(ray);
		}
	}
}


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

	AddSubdivisionStats(mBvhStats.Leaves(),
						costDecr,
						mTotalCost
						);
}


void BvHierarchy::AddSubdivisionStats(const int leaves,
									  const float renderCostDecr,
									  const float totalRenderCost)
{
	mSubdivisionStats 
			<< "#Leaves\n" << leaves << endl
			<< "#RenderCostDecrease\n" << renderCostDecr << endl 
			<< "#TotalRenderCost\n" << totalRenderCost << endl;
			//<< "#AvgRenderCost\n" << avgRenderCost << 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->mVssRays.end();

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

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


void BvHierarchy::MailViewCells(const ObjectContainer &objects, 
								const bool isFront,
								ViewCellContainer &collectedViewCells) const
{
	VssRayContainer rays;
	CollectRays(objects, rays);
	
	VssRayContainer::const_iterator rit, rit_end = rays.end();

	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);

		// view cell is assigned to left and / or right child
		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);
			
		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// mail view cell
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
				collectedViewCells.push_back(vc);

			if (isFront)
			{
				if (!vc->Mailed())
					vc->Mail();
			}
			else
			{
				if (!vc->Mailed())
					vc->Mail(1);
				else
					vc->Mail(2);
			}
		}				
	}
}


float BvHierarchy::EvalRenderCostDecrease(const ObjectContainer &objectsFront,
										  const ObjectContainer &objectsBack,
										  const BvhTraversalData &tData,
										  float &normalizedOldRenderCost) const
{
	// probability that view point lies in back / front node
	float pOverall = 0;
	float pFront = 0;
	float pBack = 0;
	float pFrontAndBack = 0;

	const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();

	BvhLeaf *leaf = tData.mNode;
	
	// sum up volume seen from the objects of left and right children
	// => the volume is the weight for the render cost equation
	ViewCell::NewMail(3);

	ViewCellContainer collectedViewCells;
	MailViewCells(objectsFront, true, collectedViewCells);
	MailViewCells(objectsBack, false, collectedViewCells);

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

	// evaluate view cells volume contribution with respect to the mail id
	for (vit = collectedViewCells.begin(); vit != vit_end; ++ vit)
	{
		AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
	}

	const int totalObjects = (int)leaf->mObjects.size();
	const int nObjectsFront = (int)objectsFront.size();
	const int nObjectsBack = (int)objectsBack.size();

	// these are non-overlapping sets
	pOverall = pFront + pBack + pFrontAndBack;

	//-- pvs rendering heuristics
	const float oldRenderCost = pOverall * totalObjects;
	const float newRenderCost = nObjectsFront * pFront + nObjectsBack * pBack + totalObjects * pFrontAndBack;

	// normalize volume with view space volume
	const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;

	Debug << "\n(((( eval render cost decrease ))))" << endl
		  << "back objects: " << nObjectsBack << " front objects " << nObjectsFront << " total objects: " << totalObjects << endl 
		  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol 
		  << " front and back p " << pFrontAndBack / viewSpaceVol << " p: " << pOverall / viewSpaceVol << endl
		  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << renderCostDecrease << endl;

	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

	if (oldRenderCost < newRenderCost * 0.99)
		Debug << "\nwarning2!!:\n" << "old rc: " << oldRenderCost * viewSpaceVol << " new rc: " << newRenderCost * viewSpaceVol << endl;
	
	return renderCostDecrease;
}


AxisAlignedBox3 BvHierarchy::ComputeBoundingBox(const ObjectContainer &objects)
{
	AxisAlignedBox3 box;
	box.Initialize();

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

	//-- compute bounding box
	for (oit = objects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

		// compute bounding box of view space
		mBoundingBox.Include(obj->GetBox());
	}
	return box;
	
}


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

	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());
		}
	}
}


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


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

	ViewCell::NewMail();

	// loop through all object and collect view cell pvs of this node
	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

		ViewCellPvsMap::const_iterator vit, vit_end = obj->mViewCellPvs.mEntries.end();

		for (vit = obj->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
		{
            ViewCell *vc = (*vit).first;

			if (!vc->Mailed())
			{
				vc->Mail();
				viewCells.push_back(vc);
			}
		}
	}
}


void BvHierarchy::CollectDirtyCandidates(BvhSubdivisionCandidate *sc, 
										 vector<SubdivisionCandidate *> &dirtyList)
{
	BvhTraversalData &tData = sc->mParentData;
	BvhLeaf *node = tData.mNode;
	
	ViewCell::NewMail();

	ViewCellContainer viewCells;
	ViewCellContainer tmpViewCells;

	VssRayContainer rays;
	CollectRays(node->mObjects, rays);

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

	// find all view cells associated with the rays, add them to view cells
	for (rit = rays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit);
		mVspTree->GetViewCells(*ray, tmpViewCells);

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

		for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
		{
			VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);

			if (!vc->Mailed())
			{
				vc->Mail();
				viewCells.push_back(vc);
			}
		}
	}


	// split candidates holding this 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 = dynamic_cast<VspViewCell *>(*vit);

		VspLeaf *leaf = vc->mLeaf;
		dirtyList.push_back(leaf->GetSubdivisionCandidate());
	}
}


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;
	}
}


BvhLeaf *BvHierarchy::GetLeaf(Intersectable *object, BvhNode *node) const
{
	// rather use the simple version
	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 = dynamic_cast<BvhLeaf *>(node);

			if (IsObjectInLeaf(leaf, object))
				return leaf;
		} 
		else 	
		{	
			// find point
			BvhInterior *interior = dynamic_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;
}


BvhIntersectable *BvHierarchy::GetOrCreateBvhIntersectable(BvhNode *node)
{
	// search nodes
	std::map<BvhNode *, BvhIntersectable *>::
		const_iterator it = mBvhIntersectables.find(node);

	if (it != mBvhIntersectables.end()) 
	{
		return (*it).second;
	}

	// not in map => create new entry
	BvhIntersectable *bvhObj = new BvhIntersectable(node);
	mBvhIntersectables[node] = bvhObj;

	return bvhObj;
}


void BvHierarchy::AddViewCellVolumeContri(ViewCell *vc,
										  float &frontVol,
										  float &backVol,
										  float &frontAndBackVol) const
{
	if (vc->Mailed())
	{
		frontVol += vc->GetVolume();
	}
	else if (vc->Mailed(1))
	{
		backVol += vc->GetVolume();
	}
	else if (vc->Mailed(2))
	{
		frontAndBackVol += vc->GetVolume();
	}
}

/*
int BvHierarchy::UpdateViewCellsPvs(BvhLeaf *leaf, 
								const RayInfoContainer &rays) const

{
	MailablePvsData::NewMail();
	
	ViewCellContainer touchedViewCells;
	CollectTouchedViewCells(rays, touchedViewCells);

	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			float contri;
			AddViewCellToObjectPvs(obj, vc, contri, true);
		}
	}

	return 0;
}


int BvHierarchy::RemoveParentViewCellsPvs(BvhLeaf *leaf, 
									  const RayInfoContainer &rays
									  ) const

{
	MailablePvsData::NewMail();
	
	ViewCellContainer touchedViewCells;
	CollectTouchedViewCells(rays, touchedViewCells);

	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

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

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			
			MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

			if (vdata && !vdata->Mailed())
			{
				vdata->Mail();
				obj->mViewCellPvs.RemoveSample(vc, 1);
			}
		}
	}

	return 0;
}
*/

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 = dynamic_cast<BvhLeaf *>(node);

		stream << "<Leaf"
			   << " min=\"" << leaf->GetBoundingBox().Min()
			   << " max=\"" << leaf->GetBoundingBox().Max() 
			   << " objects=\"";
		
		//-- export objects
		ExportObjects(leaf, stream);
		
		stream << "\" />" << endl;
	}
	else
	{	
		BvhInterior *interior = dynamic_cast<BvhInterior *>(node);
	
		stream << "<Interior"
			   << " min=\"" << interior->GetBoundingBox().Min()
			   << " max=\"" << interior->GetBoundingBox().Max()
			   << "\">" << endl;

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

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


float BvHierarchy::EvalViewCellsVolume(BvhLeaf *leaf) const
{
	float vol = 0;

	VssRayContainer rays;
	CollectRays(leaf->mObjects, rays);

	ViewCell::NewMail();

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

	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit);
	
		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);
    	
		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

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

			if (!vc->Mailed())
			{
				//Debug << "single vol: "<< vc->GetVolume() << endl;
				vc->Mail();
				vol += vc->GetVolume();
			}
		}
	}

	return vol;
}


SubdivisionCandidate *BvHierarchy::PrepareConstruction(const VssRayContainer &sampleRays,
													   ObjectContainer &objects,
													   AxisAlignedBox3 *forcedObjectSpace
													   //, RayInfoContainer &rays
													   )
{
	// store pointer to this tree
	BvhSubdivisionCandidate::sBvHierarchy = this;
	mBvhStats.nodes = 1;

	// compute bounding box from objects
	if (forcedObjectSpace)
		mBoundingBox = *forcedObjectSpace;
	else
		mBoundingBox = ComputeBoundingBox(objects);

	mTermMinVolume *= mBoundingBox.GetVolume();
	mGlobalCostMisses = 0;

	// sort objects by id in order to have sorted objects in
	// the leaf nodes
	stable_sort(objects.begin(), objects.end(), ilt);

	//-- create new root

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

	// only rays intersecting objects in node are interesting
	AssociateObjectsWithRays(sampleRays);
	// associate root with current objects
	AssociateObjectsWithLeaf(bvhleaf);

	//-- add first candidate for object space partition

	// create osp traversal data
	BvhTraversalData oData(bvhleaf, 0, mBoundingBox.GetVolume());

	//-- first split candidate
	BvhSubdivisionCandidate *oSubdivisionCandidate = 
		new BvhSubdivisionCandidate(oData);

	//UpdateViewCellsPvs(kdleaf, rays);

	EvalSubdivisionCandidate(*oSubdivisionCandidate);

// probabilty is voume of all "seen" view cells
#if 1
	const float prop = EvalViewCellsVolume(bvhleaf);
#else
	const float prop = GetBoundingBox().GetVolume();
#endif
	mTotalCost = (float)objects.size() * prop / 
		mVspTree->GetBoundingBox().GetVolume();

	EvalSubdivisionStats(*oSubdivisionCandidate);

	return oSubdivisionCandidate;
}


bool BvHierarchy::AddLeafToPvs(BvhLeaf *leaf, 
							   ViewCell *vc, 
							   const float pdf, 
							   float &contribution)
{
	// add kd intersecable to pvs
	BvhIntersectable *bvhObj = GetOrCreateBvhIntersectable(leaf);
	
	return vc->AddPvsSample(bvhObj, pdf, contribution);
}


void BvhStatistics::Print(ostream &app) const
{
	app << "=========== OspTree 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_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
		<<	maxDepthNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n" 
		<< minPvsNodes * 100 / (double)Leaves() << 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_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_INVALIDLEAVES (number of invalid leaves )\n" << invalidLeaves << endl;

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

	//app << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
	
	app << "========== END OF VspTree statistics ==========\n";
}


}