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

#include "ViewCell.h"
#include "Plane3.h"
#include "VspTree.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 "KdTree.h"
#include "IntersectableWrapper.h"
#include "HierarchyManager.h"
#include "BvHierarchy.h"
#include "OspTree.h"



namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0

/////////////
//-- static members

VspTree *VspTree::VspSubdivisionCandidate::sVspTree = NULL;
int VspNode::sMailId = 1;

// variable for debugging volume contribution for heuristics
static float debugVol;


// pvs penalty can be different from pvs size
inline static float EvalPvsPenalty(const float pvs, 
								   const float lower,
								   const float upper)
{
	// clamp to minmax values
	if (pvs < lower)
	{
		return (float)lower;
	}
	else if (pvs > upper)
	{
		return (float)upper;
	}
	return (float)pvs;
}

#if WORK_WITH_VIEWCELLS
static bool ViewCellHasMultipleReferences(Intersectable *obj, 
										  ViewCell *vc, 
										  bool checkOnlyMailed)
{
	MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

	if (vdata)
	{
		// more than one view cell sees this object inside different kd cells
		if (!checkOnlyMailed || !vdata->Mailed())
		{
			if (checkOnlyMailed)
				vdata->Mail();
			//Debug << "sumpdf: " << vdata->mSumPdf << endl;
			if (vdata->mSumPdf > 1.5f) 
				return true;
		}
	}
	
	return false;
}


void VspTree::RemoveParentViewCellReferences(ViewCell *parent) const
{
	KdLeaf::NewMail();

	// remove the parents from the object pvss
	ObjectPvsEntries::const_iterator oit, oit_end = parent->GetPvs().mEntries.end();

	for (oit = parent->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *object = (*oit).first;
		// HACK: make sure that the view cell is removed from the pvs
		const float high_contri = 9999999;

		// remove reference count of view cells
		object->mViewCellPvs.RemoveSample(parent, high_contri);
	}
}


void VspTree::AddViewCellReferences(ViewCell *vc) const
{
	KdLeaf::NewMail();

	// Add front view cell to the object pvsss
	ObjectPvsEntries::const_iterator oit, oit_end = vc->GetPvs().mEntries.end();

	for (oit = vc->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *object = (*oit).first;

		// increase reference count of view cells
		object->mViewCellPvs.AddSample(vc, 1);
	}
}

#endif

void VspTreeStatistics::Print(ostream &app) const
{
	app << "=========== VspTree 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 << "#N_SPLITS ( Number of splits in axes x y z)\n";

	for (int i = 0; i < 3; ++ i)
		app << splits[i] << " ";

	app << 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_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n" 
		<< minRaysNodes * 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_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n" 
		<<	maxRayContribNodes * 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 << "#AVGRAYS (number of rays / leaf)\n" << AvgRays() << endl;
	
	app << "#N_GLOBALCOSTMISSES ( Global cost misses )\n" << mGlobalCostMisses << endl;

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



/******************************************************************/
/*                  class VspNode implementation                  */
/******************************************************************/


VspNode::VspNode(): 
mParent(NULL), 
mTreeValid(true), 
mTimeStamp(0)
//,mRenderCostDecr(0)
//,mMemoryIncr(0)
//,mPvsEntriesIncr(0)
{}


VspNode::VspNode(VspInterior *parent): 
mParent(parent), 
mTreeValid(true), 
//,mMemoryIncr(0)
//,mRenderCostDecr(0)
//,mPvsEntriesIncr(0)
mTimeStamp(0)
{}


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


VspInterior *VspNode::GetParent() 
{ 
	return mParent; 
}


void VspNode::SetParent(VspInterior *parent)
{
	mParent = parent;
}


bool VspNode::IsSibling(VspNode *n) const
{
	return  ((this != n) && mParent && 
			 (mParent->GetFront() == n) || (mParent->GetBack() == n));
}


int VspNode::GetDepth() const
{
	int depth = 0;
	VspNode *p = mParent;
	
	while (p)
	{
		p = p->mParent;
		++ depth;
	}

	return depth;
}


bool VspNode::TreeValid() const
{
	return mTreeValid;
}


void VspNode::SetTreeValid(const bool v)
{
	mTreeValid = v;
}



/****************************************************************/
/*              class VspInterior implementation                */
/****************************************************************/


VspInterior::VspInterior(const AxisAlignedPlane &plane): 
mPlane(plane), mFront(NULL), mBack(NULL)
{}


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


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


VspNode *VspInterior::GetBack() 
{
	return mBack;
}


VspNode *VspInterior::GetFront() 
{
	return mFront;
}


AxisAlignedPlane VspInterior::GetPlane() const
{
	return mPlane;
}


float VspInterior::GetPosition() const
{
	return mPlane.mPosition;
}


int VspInterior::GetAxis() const
{
	return mPlane.mAxis;
}


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


void VspInterior::SetupChildLinks(VspNode *front, VspNode *back) 
{
    mBack = back;
    mFront = front;
}


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


void VspInterior::SetBoundingBox(const AxisAlignedBox3 &box)
{
	mBoundingBox = box;
}


int VspInterior::Type() const
{
	return Interior;
}



/****************************************************************/
/*                  class VspLeaf implementation                */
/****************************************************************/


VspLeaf::VspLeaf(): 
mViewCell(NULL), mSubdivisionCandidate(NULL)
{
}


VspLeaf::~VspLeaf()
{
	VssRayContainer::const_iterator vit, vit_end = mVssRays.end();
	for (vit = mVssRays.begin(); vit != vit_end; ++ vit)
	{
		VssRay *ray = *vit;
		ray->Unref();

		if (!ray->IsActive())
			delete ray;
	}
	//CLEAR_CONTAINER(mVssRays);
}


int VspLeaf::Type() const
{
	return Leaf;
}


VspLeaf::VspLeaf(ViewCellLeaf *viewCell): 
mViewCell(viewCell)
{
}


VspLeaf::VspLeaf(VspInterior *parent): 
VspNode(parent), mViewCell(NULL)
{}


VspLeaf::VspLeaf(VspInterior *parent, ViewCellLeaf *viewCell): 
VspNode(parent), mViewCell(viewCell)
{
}


ViewCellLeaf *VspLeaf::GetViewCell() const
{
	return mViewCell;
}


void VspLeaf::SetViewCell(ViewCellLeaf *viewCell)
{
	mViewCell = viewCell;
}


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



/*************************************************************************/
/*                       class VspTree implementation                    */
/*************************************************************************/


VspTree::VspTree():
mRoot(NULL),
mOutOfBoundsCell(NULL),
mStoreRays(false),
mTimeStamp(1),
mHierarchyManager(NULL)
{
	mLocalSubdivisionCandidates = new vector<SortableEntry>;

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

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

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

	Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("VspTree.Termination.minPvs", mTermMinPvs);
	Environment::GetSingleton()->GetIntValue("VspTree.Termination.minRays", mTermMinRays);
	Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minProbability", mTermMinProbability);
	Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxRayContribution", mTermMaxRayContribution);
	
	Environment::GetSingleton()->GetIntValue("VspTree.Termination.missTolerance", mTermMissTolerance);
	Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxViewCells", mMaxViewCells);
	// max cost ratio for early tree termination
	Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxCostRatio", mTermMaxCostRatio);

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

	Environment::GetSingleton()->GetFloatValue("VspTree.maxStaticMemory", mMaxMemory);


	//////////////
	//-- factors for bsp tree split plane heuristics

	Environment::GetSingleton()->GetFloatValue("VspTree.Termination.ct_div_ci", mCtDivCi);
	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.epsilon", mEpsilon);
	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.minBand", mMinBand);
	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.maxBand", mMaxBand);
	Environment::GetSingleton()->GetIntValue("VspTree.maxTests", mMaxTests);

	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
	
	// if only the driving axis is used for axis aligned split
	Environment::GetSingleton()->GetBoolValue("VspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
	Environment::GetSingleton()->GetBoolValue("VspTree.useCostHeuristics", mUseCostHeuristics);
	Environment::GetSingleton()->GetBoolValue("VspTree.simulateOctree", mCirculatingAxis);

	//mMemoryConst = (float)(sizeof(VspLeaf) + sizeof(VspViewCell));
	//mMemoryConst = 50;//(float)(sizeof(VspViewCell));

	//mMemoryConst = (float)(sizeof(VspLeaf) + sizeof(ObjectPvs));
	mMemoryConst = 16;//(float)sizeof(ObjectPvs);
	cout << "vsp memcost: " << mMemoryConst << endl;

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

	Debug << "******* VSP options ******** " << endl;

    Debug << "max depth: " << mTermMaxDepth << endl;
	Debug << "min PVS: " << mTermMinPvs << endl;
	Debug << "min probabiliy: " << mTermMinProbability << endl;
	Debug << "min rays: " << mTermMinRays << endl;
	Debug << "max ray contri: " << mTermMaxRayContribution << endl;
	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	Debug << "miss tolerance: " << mTermMissTolerance << endl;
	Debug << "max view cells: " << mMaxViewCells << 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 << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;

	Debug << "circulating axis: " << mCirculatingAxis << endl;
	Debug << "minband: " << mMinBand << endl;
	Debug << "maxband: " << mMaxBand << endl;

	Debug << "vsp mem const: " << mMemoryConst << endl;

	cout << "use cost heuristics: " << mUseCostHeuristics << endl;

	Debug << endl;
}


VspViewCell *VspTree::GetOutOfBoundsCell()
{
	return mOutOfBoundsCell;
}


VspViewCell *VspTree::GetOrCreateOutOfBoundsCell()
{
	if (!mOutOfBoundsCell)
	{
		mOutOfBoundsCell = new VspViewCell();
		mOutOfBoundsCell->SetId(-1);
		mOutOfBoundsCell->SetValid(false);
	}

	return mOutOfBoundsCell;
}


const VspTreeStatistics &VspTree::GetStatistics() const
{
	return mVspStats;
}


VspTree::~VspTree()
{
	DEL_PTR(mRoot);
	DEL_PTR(mLocalSubdivisionCandidates);
}


void VspTree::ComputeBoundingBox(const VssRayContainer &rays,
								 AxisAlignedBox3 *forcedBoundingBox) 
{	
	if (forcedBoundingBox)
	{
		mBoundingBox = *forcedBoundingBox;
		return;
	}
	
	//////////////////////////////////////////////
	// bounding box of view space includes all visibility events
	mBoundingBox.Initialize();
	VssRayContainer::const_iterator rit, rit_end = rays.end();

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

		mBoundingBox.Include(ray->GetTermination());
		mBoundingBox.Include(ray->GetOrigin());
	}
}


void VspTree::AddSubdivisionStats(const int viewCells,
								  const float renderCostDecr,
								  const float totalRenderCost,
								  const float avgRenderCost)
{
	mSubdivisionStats 
			<< "#ViewCells\n" << viewCells << endl
			<< "#RenderCostDecrease\n" << renderCostDecr << endl 
			<< "#TotalRenderCost\n" << totalRenderCost << endl
			<< "#AvgRenderCost\n" << avgRenderCost << endl;
}


// $$matt temporary: number of rayrefs + pvs should be in there, but
// commented out for testing
float VspTree::GetMemUsage() const
{
	return (float)
		 (sizeof(VspTree)
		  + mVspStats.Leaves() * sizeof(VspLeaf) 
		  + mVspStats.Leaves() * sizeof(VspViewCell)
		  + mVspStats.Interior() * sizeof(VspInterior)
		  //+ mVspStats.pvs * sizeof(PvsData)
		  //+ mVspStats.rayRefs * sizeof(RayInfo)
		  ) / (1024.0f * 1024.0f);
}


inline bool VspTree::LocalTerminationCriteriaMet(const VspTraversalData &data) const
{
	const bool localTerminationCriteriaMet = (0
		|| ((int)data.mRays->size() <= mTermMinRays)
		|| (data.mPvs <= mTermMinPvs)
		|| (data.mProbability <= mTermMinProbability)
		//|| (data.GetAvgRayContribution() > mTermMaxRayContribution)
		|| (data.mDepth >= mTermMaxDepth)
		);

#if GTP_DEBUG
	if (localTerminationCriteriaMet)
	{
		Debug << "local termination criteria met:" << endl;
		Debug << "rays: " << (int)data.mRays->size() << "  " << mTermMinRays << endl;
		Debug << "pvs: " << data.mPvs << " " << mTermMinPvs << endl;
		Debug << "p: " <<  data.mProbability << " " << mTermMinProbability << endl;
		Debug << "avg contri: " << data.GetAvgRayContribution() << " " << mTermMaxRayContribution << endl;
		Debug << "depth " << data.mDepth << " " << mTermMaxDepth << endl;
	}
#endif
	return localTerminationCriteriaMet;		
}


inline bool VspTree::GlobalTerminationCriteriaMet(const VspTraversalData &data) const
{
	// note: to track for global cost misses does not really 
	// make sense because cost termination  happens in the hierarchy mananger

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

#if GTP_DEBUG
	if (terminationCriteriaMet)
	{
		Debug << "vsp global termination criteria met:" << endl;
		Debug << "cost misses: " << mVspStats.mGlobalCostMisses << " " << mTermGlobalCostMissTolerance << endl;
		Debug << "leaves: " << mVspStats.Leaves() << " " <<  mMaxViewCells << endl;
	}
#endif

	return terminationCriteriaMet;
}


void VspTree::CreateViewCell(VspTraversalData &tData, const bool updatePvs)
{
	///////////////
	//-- create new view cell

	VspLeaf *leaf = tData.mNode;

	VspViewCell *viewCell = new VspViewCell();
    leaf->SetViewCell(viewCell);
	
	int conSamp = 0;
	float sampCon = 0.0f;

	if (updatePvs)
	{
		// update pvs of view cell
		AddSamplesToPvs(leaf, *tData.mRays, sampCon, conSamp);

		// update scalar pvs size value
		ObjectPvs &pvs = viewCell->GetPvs();
		mViewCellsManager->UpdateScalarPvsSize(viewCell, pvs.EvalPvsCost(), pvs.GetSize());

		mVspStats.contributingSamples += conSamp;
		mVspStats.sampleContributions += (int)sampCon;
	}

	if (mStoreRays)
	{
		///////////
		//-- store sampling rays

        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();

		for (it = tData.mRays->begin(); it != it_end; ++ it)
		{
			(*it).mRay->Ref();			
			// note: should rather store rays with view cell
			leaf->mVssRays.push_back((*it).mRay);
		}
	}
		
	// set view cell values
	viewCell->mLeaves.push_back(leaf);

	viewCell->SetVolume(tData.mProbability);
    leaf->mProbability = tData.mProbability;
}


void VspTree::EvalSubdivisionStats(const SubdivisionCandidate &sc)
{
	const float costDecr = sc.GetRenderCostDecrease();
	
	AddSubdivisionStats(mVspStats.Leaves(),
						costDecr,
						mTotalCost,
						(float)mTotalPvsSize / (float)mVspStats.Leaves());
}


VspNode *VspTree::Subdivide(SplitQueue &tQueue,
							SubdivisionCandidate *splitCandidate,
							const bool globalCriteriaMet)
{
	// todo remove dynamic cast
	VspSubdivisionCandidate *sc = 
		dynamic_cast<VspSubdivisionCandidate *>(splitCandidate);

	VspTraversalData &tData = sc->mParentData;
	VspNode *newNode = tData.mNode;

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

		VspTraversalData tFrontData;
		VspTraversalData tBackData;
		
		// create new interior node and two leaf node
		const int maxCostMisses = sc->GetMaxCostMisses();

		newNode = SubdivideNode(*sc,
								tFrontData, 
								tBackData);
	
		// how often was max cost ratio missed in this branch?
		tFrontData.mMaxCostMisses = maxCostMisses;
		tBackData.mMaxCostMisses = maxCostMisses;
			
		mTotalCost -= sc->GetRenderCostDecrease();
		//mTotalPvsSize += (int)(tFrontData.mPvs + tBackData.mPvs - tData.mPvs);
		mPvsEntries += sc->GetPvsEntriesIncr();

		// subdivision statistics
		if (1) EvalSubdivisionStats(*sc);
		
		/////////////
		//-- evaluate new split candidates for global greedy cost heuristics

		VspSubdivisionCandidate *frontCandidate = new VspSubdivisionCandidate(tFrontData);
		VspSubdivisionCandidate *backCandidate = new VspSubdivisionCandidate(tBackData);

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

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

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

		// note: leaf is not destroyed because it is needed to collect 
		// dirty candidates in hierarchy manager
	}

	if (newNode->IsLeaf()) // subdivision terminated
	{
		VspLeaf *leaf = dynamic_cast<VspLeaf *>(newNode);
		
#if 0
        /////////////
		//-- store pvs optained from rays

		// view cell is created during subdivision
		ViewCell *viewCell = leaf->GetViewCell();

		int conSamp = 0;
		float sampCon = 0.0f;

		AddSamplesToPvs(leaf, *tData.mRays, sampCon, conSamp);

		// update scalar pvs size value
		ObjectPvs &pvs = viewCell->GetPvs();
		mViewCellsManager->UpdateScalarPvsSize(viewCell, pvs.EvalPvsCost(), pvs.GetSize());

		mVspStats.contributingSamples += conSamp;
		mVspStats.sampleContributions += (int)sampCon;
#endif
		if (mStoreRays)
		{
			//////////
			//-- store rays piercing this view cell
			RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
			for (it = tData.mRays->begin(); it != it_end; ++ it)
			{
				(*it).mRay->Ref();			
				leaf->mVssRays.push_back((*it).mRay);
				//leaf->mVssRays.push_back(new VssRay(*(*it).mRay));
			}
		}

		// finally evaluate statistics for this leaf
		EvaluateLeafStats(tData);
		// detach subdivision candidate: this leaf is no candidate for
		// splitting anymore
		tData.mNode->SetSubdivisionCandidate(NULL); 
		// detach node so it won't get deleted
		tData.mNode = NULL;
	}

	return newNode;
}


void VspTree::EvalSubdivisionCandidate(VspSubdivisionCandidate &splitCandidate, 
									   bool computeSplitPlane)
{
	if (computeSplitPlane)
	{
		float frontProb;
		float backProb;

		// compute locally best split plane
		const float ratio = SelectSplitPlane(splitCandidate.mParentData, 
									 		 splitCandidate.mSplitPlane, 
											 frontProb, 
											 backProb);
	
		const bool maxCostRatioViolated = mTermMaxCostRatio < ratio;

		const int maxCostMisses = splitCandidate.mParentData.mMaxCostMisses;
		// max cost threshold violated?
		splitCandidate.SetMaxCostMisses(maxCostRatioViolated ? 
										maxCostMisses + 1: maxCostMisses);
	}
	
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(splitCandidate.mParentData.mNode);

	/////////////
	// avg ray contri

	const int pvs = EvalPvsEntriesSize(*splitCandidate.mParentData.mRays);

	const float avgRayContri = (float)pvs / 
		((float)splitCandidate.mParentData.mRays->size() + Limits::Small);

	splitCandidate.SetAvgRayContribution(avgRayContri);
	
	// compute global decrease in render cost
	float oldRenderCost;
	const float renderCostDecr = EvalRenderCostDecrease(splitCandidate, oldRenderCost);
	splitCandidate.SetRenderCostDecrease(renderCostDecr);

	// the increase in pvs entries num induced by this split
	const int pvsEntriesIncr = EvalPvsEntriesIncr(splitCandidate);
	splitCandidate.SetPvsEntriesIncr(pvsEntriesIncr);

	// take render cost of node into account 
	// otherwise danger of being stuck in a local minimum!
	const float factor = mRenderCostDecreaseWeight;

	float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;

	if (mHierarchyManager->mConsiderMemory)
	{
		priority /= ((float)splitCandidate.GetPvsEntriesIncr() + mMemoryConst);
	}

	splitCandidate.SetPriority(priority);
}


int VspTree::EvalPvsEntriesIncr(VspSubdivisionCandidate &splitCandidate) const
{
	float oldPvsSize = 0;
	float fPvsSize = 0;
	float bPvsSize = 0;
	
	const AxisAlignedPlane candidatePlane = splitCandidate.mSplitPlane;
	
	Intersectable::NewMail(3);
	KdLeaf::NewMail(3);

	RayInfoContainer::const_iterator rit, rit_end = splitCandidate.mParentData.mRays->end();

    // this is the main ray classification loop!
	for(rit = splitCandidate.mParentData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
		RayInfo rayInf = *rit;
		
		float t;
		// classify ray
		const int cf = 	rayInf.ComputeRayIntersection(candidatePlane.mAxis, 
													  candidatePlane.mPosition, t);

		UpdatePvsEntriesContribution(*ray, true, cf, fPvsSize, bPvsSize, oldPvsSize);
#if COUNT_ORIGIN_OBJECTS
		UpdatePvsEntriesContribution(*ray, false, cf, fPvsSize, bPvsSize, oldPvsSize);
#endif
	}

	const float oldPvsRatio = (splitCandidate.mParentData.mPvs > 0) ? 
		oldPvsSize / splitCandidate.mParentData.mPvs : 1;
	const float correctedOldPvs = splitCandidate.mParentData.mCorrectedPvs * oldPvsRatio;

	splitCandidate.mCorrectedFrontPvs = 
		mHierarchyManager->EvalCorrectedPvs(fPvsSize, 
											correctedOldPvs, 
											splitCandidate.GetAvgRayContribution());
	splitCandidate.mCorrectedBackPvs = 
		mHierarchyManager->EvalCorrectedPvs(bPvsSize, 
											correctedOldPvs, 
											splitCandidate.GetAvgRayContribution());
	
	if (0)
	cout << "vsp pvs"
		 << " avg ray contri: " << splitCandidate.GetAvgRayContribution() << " ratio: " << oldPvsRatio 
		 << " parent: " << correctedOldPvs << " " << " old vol: " << oldPvsSize
		 << " frontpvs: " << fPvsSize << " corr. " << splitCandidate.mCorrectedFrontPvs 
		 << " backpvs: " << bPvsSize << " corr. " << splitCandidate.mCorrectedBackPvs << endl;


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


VspInterior *VspTree::SubdivideNode(const VspSubdivisionCandidate &sc,
									VspTraversalData &frontData,
									VspTraversalData &backData)
{
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(sc.mParentData.mNode);

	const VspTraversalData &tData = sc.mParentData;
	const AxisAlignedPlane splitPlane = sc.mSplitPlane;

	///////////////
	//-- new traversal values

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

	frontData.mRays = new RayInfoContainer();
	backData.mRays = new RayInfoContainer();

	//-- subdivide rays
	SplitRays(splitPlane, *tData.mRays, *frontData.mRays, *backData.mRays);

	//-- compute pvs
	frontData.mPvs = (float)EvalPvsEntriesSize(*frontData.mRays);
	backData.mPvs = (float)EvalPvsEntriesSize(*backData.mRays);
	
	//-- compute pvs correction for coping with undersampling
	frontData.mCorrectedPvs = sc.mCorrectedFrontPvs;
	backData.mCorrectedPvs = sc.mCorrectedBackPvs;

	//-- split front and back node geometry and compute area
	tData.mBoundingBox.Split(splitPlane.mAxis, 
							 splitPlane.mPosition, 
							 frontData.mBoundingBox, 
							 backData.mBoundingBox);

	frontData.mProbability = frontData.mBoundingBox.GetVolume();
	backData.mProbability = tData.mProbability - frontData.mProbability;

	//-- compute render cost
	frontData.mRenderCost = (float)EvalPvsCost(*frontData.mRays);
	backData.mRenderCost = (float)EvalPvsCost(*backData.mRays);
	
	frontData.mCorrectedRenderCost = sc.mCorrectedFrontRenderCost;
	backData.mCorrectedRenderCost = sc.mCorrectedBackRenderCost;

	////////
	//-- update some stats
	
	if (tData.mDepth > mVspStats.maxDepth)
	{	
		mVspStats.maxDepth = tData.mDepth;
	}

	// two more leaves per split
	mVspStats.nodes += 2;
	// and a new split
	++ mVspStats.splits[splitPlane.mAxis];


	////////////////
	//-- create front and back and subdivide further

	VspInterior *interior = new VspInterior(splitPlane);
	VspInterior *parent = leaf->GetParent();

	// replace a link from node's parent
	if (parent)
	{
		parent->ReplaceChildLink(leaf, interior);
		interior->SetParent(parent);

#if WORK_WITH_VIEWCELLS
		// remove "parent" view cell from pvs of all objects (traverse through rays)
		RemoveParentViewCellReferences(tData.mNode->GetViewCell());
#endif
	}
	else // new root
	{
		mRoot = interior;
	}

	VspLeaf *frontLeaf = new VspLeaf(interior);
	VspLeaf *backLeaf = new VspLeaf(interior);

	// and setup child links
	interior->SetupChildLinks(frontLeaf, backLeaf);
	
	// add bounding box
	interior->SetBoundingBox(tData.mBoundingBox);

	// set front and back leaf
	frontData.mNode = frontLeaf;
	backData.mNode = backLeaf;

	// create front and back view cell for the new leaves
	CreateViewCell(frontData, false);
	CreateViewCell(backData, false);

	// set the time stamp so the order of traversal can be reconstructed
	interior->mTimeStamp = mHierarchyManager->mTimeStamp ++;

#if WORK_WITH_VIEWCELL_PVS
	// create front and back view cell
	// add front and back view cell to 
	// "potentially visible view cells" 
	// of the objects in front and back pvs

	AddViewCellReferences(frontLeaf->GetViewCell());
	AddViewCellReferences(backLeaf->GetViewCell());
#endif

	return interior;
}


void VspTree::AddSamplesToPvs(VspLeaf *leaf,
							  const RayInfoContainer &rays,
							  float &sampleContributions,
							  int &contributingSamples)
{
	sampleContributions = 0;
	contributingSamples = 0;
  
	RayInfoContainer::const_iterator it, it_end = rays.end();
  
	ViewCellLeaf *vc = leaf->GetViewCell();

	// add contributions from samples to the pvs
	for (it = rays.begin(); it != it_end; ++ it)
	{
		float sc = 0.0f;
		VssRay *ray = (*it).mRay;

		bool madeContrib = false;
		float contribution = 1;

		Intersectable *entry = 
			mHierarchyManager->GetIntersectable(ray->mTerminationObject, true);

		if (entry) 
		{
			madeContrib = 
				vc->GetPvs().AddSample(entry, ray->mPdf);

			sc += contribution;
		}
#if COUNT_ORIGIN_OBJECTS
		entry = mHierarchyManager->GetIntersectable(ray->mOriginObject, true);

		if (entry) 
		{
			madeContrib = 
				vc->GetPvs().AddSample(entry, ray->mPdf);

			sc += contribution;
		}
#endif
		if (madeContrib)
		{
			++ contributingSamples;
		}

		// store rays for visualization
		if (0) leaf->mVssRays.push_back(new VssRay(*ray));
	}
}


void VspTree::SortSubdivisionCandidates(const RayInfoContainer &rays,
								  const int axis, 
								  float minBand, 
								  float maxBand)
{
	mLocalSubdivisionCandidates->clear();

	const int requestedSize = 2 * (int)(rays.size());

	// creates a sorted split candidates array
	if (mLocalSubdivisionCandidates->capacity() > 500000 &&
		requestedSize < (int)(mLocalSubdivisionCandidates->capacity() / 10) )
	{
    	delete mLocalSubdivisionCandidates;
   		mLocalSubdivisionCandidates = new vector<SortableEntry>;
	}

	mLocalSubdivisionCandidates->reserve(requestedSize);

	float pos;
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	const bool delayMinEvent = false;

	////////////
	//-- insert all queries
	for (rit = rays.begin(); rit != rit_end; ++ rit)
	{
		const bool positive = (*rit).mRay->HasPosDir(axis);
		
		// origin point
		pos = (*rit).ExtrapOrigin(axis);
		const int oType = positive ? SortableEntry::ERayMin : SortableEntry::ERayMax;

		if (delayMinEvent && oType == SortableEntry::ERayMin)
			pos += mEpsilon; // could be useful feature for walls

		mLocalSubdivisionCandidates->push_back(SortableEntry(oType, pos, (*rit).mRay));

		// termination point
		pos = (*rit).ExtrapTermination(axis);
		const int tType = positive ? SortableEntry::ERayMax : SortableEntry::ERayMin;

		if (delayMinEvent && tType == SortableEntry::ERayMin)
			pos += mEpsilon; // could be useful feature for walls

		mLocalSubdivisionCandidates->push_back(SortableEntry(tType, pos, (*rit).mRay));
	}

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


float VspTree::PrepareHeuristics(KdLeaf *leaf)
{	
	float pvsSize = 0;
	
	if (!leaf->Mailed())
	{
		leaf->Mail();
		leaf->mCounter = 1;
		// add objects without the objects which are in several kd leaves
		pvsSize += (int)(leaf->mObjects.size() - leaf->mMultipleObjects.size());
	}
	else
	{
		++ leaf->mCounter;
	}

	//-- the objects belonging to several leaves must be handled seperately
	ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();

	for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *object = *oit;
						
		if (!object->Mailed())
		{
			object->Mail();
			object->mCounter = 1;

			++ pvsSize;
		}
		else
		{
			++ object->mCounter;
		}
	}
	
	return pvsSize;
}


float VspTree::PrepareHeuristics(const RayInfoContainer &rays)
{	
	Intersectable::NewMail();
	KdNode::NewMail();

	float pvsSize = 0;

	RayInfoContainer::const_iterator ri, ri_end = rays.end();

    // set all kd nodes / objects as belonging to the front pvs
	for (ri = rays.begin(); ri != ri_end; ++ ri)
	{
		VssRay *ray = (*ri).mRay;
		
		pvsSize += PrepareHeuristics(*ray, true);
#if COUNT_ORIGIN_OBJECTS
		pvsSize += PrepareHeuristics(*ray, false);
#endif
	}

	return pvsSize;
}


int VspTree::EvalMaxEventContribution(KdLeaf *leaf) const
{
	int pvs = 0;

	// leaf falls out of right pvs
	if (-- leaf->mCounter == 0)
	{
		pvs -= ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
	}

	//////
	//-- separately handle objects which are in several kd leaves

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

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

		if (-- object->mCounter == 0)
		{
			++ pvs;
		}
	}

	return pvs;
}


int VspTree::EvalMinEventContribution(KdLeaf *leaf) const
{
	if (leaf->Mailed())
		return 0;
	
	leaf->Mail();

	// add objects without those which are part of several kd leaves
	int pvs = ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());

	// separately handle objects which are part of several kd leaves 
	ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();

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

		// object not previously in pvs
		if (!object->Mailed())
		{
			object->Mail();
			++ pvs;
		}
	}	

	return pvs;
}


void VspTree::EvalHeuristics(const SortableEntry &ci,
							 float &pvsLeft,
							 float &pvsRight) const
{
	VssRay *ray = ci.ray;

	// eval changes in pvs causes by min event 
	if (ci.type == SortableEntry::ERayMin)
	{
		pvsLeft += EvalMinEventContribution(*ray, true);
#if COUNT_ORIGIN_OBJECTS
		pvsLeft += EvalMinEventContribution(*ray, false);
#endif
	}
	else // eval changes in pvs causes by max event 
	{
		pvsRight -= EvalMaxEventContribution(*ray, true);
#if COUNT_ORIGIN_OBJECTS
		pvsRight -= EvalMaxEventContribution(*ray, false);
#endif
	}
}


float VspTree::EvalLocalCostHeuristics(const VspTraversalData &tData,
									   const AxisAlignedBox3 &box,
									   const int axis,
									   float &position,
									   const RayInfoContainer &usedRays)
{
	const float minBox = box.Min(axis);
	const float maxBox = box.Max(axis);

	const float sizeBox = maxBox - minBox;

	const float minBand = minBox + mMinBand * sizeBox;
	const float maxBand = minBox + mMaxBand * sizeBox;

	SortSubdivisionCandidates(usedRays, axis, minBand, maxBand);

	// prepare the sweep
	// note: returns pvs size => no need t give pvs size as function parameter
	const float pvsCost = PrepareHeuristics(usedRays);

	// go through the lists, count the number of objects left and right
	// and evaluate the following cost funcion:
	// C = ct_div_ci  + (ql*rl + qr*rr)/queries

	float pvsl = 0;
	float pvsr = pvsCost;

	float pvsBack = pvsl;
	float pvsFront = pvsr;

	float sum = pvsCost * sizeBox;
	float minSum = 1e20f;

	// if no good split can be found, take mid split
	position = minBox + 0.5f * sizeBox;
	
	// the relative cost ratio
	float ratio = 99999999.0f;
	bool splitPlaneFound = false;

	Intersectable::NewMail();
	KdLeaf::NewMail();

	const float volRatio = 
		tData.mBoundingBox.GetVolume() / (sizeBox * mBoundingBox.GetVolume());

	////////
	//-- iterate through visibility events

	vector<SortableEntry>::const_iterator ci, 
		ci_end = mLocalSubdivisionCandidates->end();

#ifdef GTP_DEBUG
	const int leaves = mVspStats.Leaves();
	const bool printStats = ((axis == 0) && (leaves > 0) && (leaves < 90));
	
	ofstream sumStats;
	ofstream pvslStats;
	ofstream pvsrStats;

	if (printStats)
	{
		char str[64]; 
		
		sprintf(str, "tmp/vsp_heur_sum-%04d.log", leaves);
		sumStats.open(str);
		sprintf(str, "tmp/vsp_heur_pvsl-%04d.log", leaves);
		pvslStats.open(str);
		sprintf(str, "tmp/vsp_heur_pvsr-%04d.log", leaves);
		pvsrStats.open(str);
	}

#endif
	for (ci = mLocalSubdivisionCandidates->begin(); ci != ci_end; ++ ci)
	{
		// compute changes to front and back pvs
		EvalHeuristics(*ci, pvsl, pvsr);

		// Note: sufficient to compare size of bounding boxes of front and back side?
		if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
		{
			float currentPos;
			
			// HACK: current positition is BETWEEN visibility events
			if (0 && ((ci + 1) != ci_end))
				currentPos = ((*ci).value + (*(ci + 1)).value) * 0.5f;
			else
				currentPos = (*ci).value;			

			sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
			
#ifdef GTP_DEBUG
			if (printStats)
			{
				sumStats
					<< "#Position\n" << currentPos << endl
					<< "#Sum\n" << sum * volRatio << endl
					<< "#Pvs\n" << pvsl + pvsr << endl;

				pvslStats
					<< "#Position\n" << currentPos << endl
					<< "#Pvsl\n" << pvsl << endl;

				pvsrStats
					<< "#Position\n" << currentPos << endl
					<< "#Pvsr\n" << pvsr << endl;
			}
#endif

			if (sum < minSum)
			{
				splitPlaneFound = true;

				minSum = sum;
				position = currentPos;
				
				pvsBack = pvsl;
				pvsFront = pvsr;
			}
		}
	}
	
	/////////	
	//-- compute cost

	const float pOverall = sizeBox;
	const float pBack = position - minBox;
	const float pFront = maxBox - position;

	const float penaltyOld = pvsCost;
    const float penaltyFront = pvsFront;
	const float penaltyBack = pvsBack;
	
	const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

	if (splitPlaneFound)
	{
		ratio = newRenderCost / oldRenderCost;
	}
	
#ifdef GTP_DEBUG
	cout << "\n((((( eval local cost )))))" << endl
		  << "back pvs: " << penaltyBack << " front pvs: " << penaltyFront << " total pvs: " << penaltyOld << endl 
		  << "back p: " << pBack * volRatio << " front p " << pFront * volRatio << " p: " << pOverall * volRatio << endl
		  << "old rc: " << oldRenderCost * volRatio << " new rc: " << newRenderCost * volRatio << endl
		  << "render cost decrease: " << oldRenderCost * volRatio - newRenderCost * volRatio << endl;
#endif
	return ratio;
}


float VspTree::SelectSplitPlane(const VspTraversalData &tData,
								AxisAlignedPlane &plane,
								float &pFront,
								float &pBack)
{
	float nPosition[3];
	float nCostRatio[3];
	float nProbFront[3];
	float nProbBack[3];

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

	// do we use some kind of specialised "fixed" axis?
    const bool useSpecialAxis = 
		mOnlyDrivingAxis || mCirculatingAxis;
	
	// get subset of rays
	RayInfoContainer randomRays;
	randomRays.reserve(min(mMaxTests, (int)tData.mRays->size()));

	RayInfoContainer *usedRays;

	if (mMaxTests < (int)tData.mRays->size())
	{
		GetRayInfoSets(*tData.mRays, mMaxTests, randomRays);
		usedRays = &randomRays;
	}
	else
	{
		usedRays = tData.mRays;
	}

	if (mCirculatingAxis)
	{
		int parentAxis = 0;
		VspNode *parent = tData.mNode->GetParent();

		if (parent)
		{
			parentAxis = dynamic_cast<VspInterior *>(parent)->GetAxis();
		}

		sAxis = (parentAxis + 1) % 3;
	}
	else if (mOnlyDrivingAxis)
	{
		sAxis = box.Size().DrivingAxis();
	}
	
	for (int axis = 0; axis < 3; ++ axis)
	{
		if (!useSpecialAxis || (axis == sAxis))
		{
			if (mUseCostHeuristics)
			{
				//-- place split plane using heuristics
				nCostRatio[axis] =
					EvalLocalCostHeuristics(tData,
											box,
											axis,
											nPosition[axis],
											*usedRays);			
			}
			else
			{
				//-- split plane position is spatial median				
				nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
				nCostRatio[axis] = EvalLocalSplitCost(tData,
													  box,
													  axis,
													  nPosition[axis],
													  nProbFront[axis], 
													  nProbBack[axis]);
			}
						
			if (bestAxis == -1)
			{
				bestAxis = axis;
			}
			else if (nCostRatio[axis] < nCostRatio[bestAxis])
			{
				bestAxis = axis;
			} 
		}
	}

	/////////////////////////
	//-- assign values of best split

	plane.mAxis = bestAxis;
	// best split plane position
	plane.mPosition = nPosition[bestAxis];

	pFront = nProbFront[bestAxis];
	pBack = nProbBack[bestAxis];

	return nCostRatio[bestAxis];
}


float VspTree::EvalRenderCostDecrease(VspSubdivisionCandidate &sc,
									  float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;
	float totalPvs = 0;

	const float viewSpaceVol = mBoundingBox.GetVolume();
	const VspTraversalData &tData = sc.mParentData;
	
	const AxisAlignedPlane &candidatePlane = sc.mSplitPlane;
	const float avgRayContri = sc.GetAvgRayContribution();

	//////////////////////////////////////////////
	// mark objects in the front / back / both using mailboxing
	// then count pvs sizes

	Intersectable::NewMail(3);
	KdLeaf::NewMail(3);

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();

	for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		RayInfo rayInf = *rit;

		float t;
		
		// classify ray
		const int cf = rayInf.ComputeRayIntersection(candidatePlane.mAxis, 
													 candidatePlane.mPosition, 
													 t);

		VssRay *ray = rayInf.mRay;

		// evaluate contribution of ray endpoint to front 
		// and back pvs with respect to the classification
		UpdateContributionsToPvs(*ray, true, cf, pvsFront, pvsBack, totalPvs);
#if COUNT_ORIGIN_OBJECTS
		UpdateContributionsToPvs(*ray, false, cf, pvsFront, pvsBack, totalPvs);
#endif
	}
    
	AxisAlignedBox3 frontBox;
	AxisAlignedBox3 backBox;

	tData.mBoundingBox.Split(candidatePlane.mAxis, 
							 candidatePlane.mPosition, 
							 frontBox, 
							 backBox);

    // probability that view point lies in back / front node
	float pOverall = tData.mProbability;
	float pFront = frontBox.GetVolume();
	float pBack = pOverall - pFront;


	///////////////////
	//-- evaluate render cost heuristics

	const float oldRenderCostRatio = (tData.mRenderCost > 0)? 
		(totalPvs / tData.mRenderCost) : 1;

	const float penaltyOld = tData.mCorrectedRenderCost * oldRenderCostRatio;

	sc.mCorrectedFrontRenderCost = mHierarchyManager->EvalCorrectedPvs(pvsFront, penaltyOld, avgRayContri);
	sc.mCorrectedBackRenderCost = mHierarchyManager->EvalCorrectedPvs(pvsBack, penaltyOld, avgRayContri);

	const float oldRenderCost = pOverall * penaltyOld;
	const float newRenderCost = sc.mCorrectedFrontRenderCost * pFront + 
								sc.mCorrectedFrontRenderCost * pBack;

	// we also return the old render cost
	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

	// the render cost decrase for this split
	const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;

	if (0)
	cout << "vsp render cost"
		 << " avg ray contri: " << avgRayContri << " ratio: " << oldRenderCostRatio
		 << " parent: " << penaltyOld << " " << " old vol: " << totalPvs 
		 << " front cost: " << sc.mCorrectedFrontRenderCost << " corr. " << sc.mCorrectedFrontRenderCost
		 << " back cost: " << sc.mCorrectedBackRenderCost << " corr. " << sc.mCorrectedBackRenderCost << endl;

#ifdef GTP_DEBUG
	Debug << "\nvsp render cost decrease" << endl
		  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl 
		  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol << " p: " << pOverall / viewSpaceVol << endl
		  << "old rc: " << normalizedOldRenderCost << " new rc: " << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << renderCostDecrease << endl;
#endif

	return renderCostDecrease;
}


float VspTree::EvalLocalSplitCost(const VspTraversalData &tData,
								  const AxisAlignedBox3 &box,
								  const int axis,
								  const float &position,
								  float &pFront,
								  float &pBack) const
{
	float pvsTotal = 0;
	float pvsFront = 0;
	float pvsBack = 0;
	
	// create unique ids for pvs heuristics
	Intersectable::NewMail(3);
	KdLeaf::NewMail(3);

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();

	// this is the main ray classification loop!
	for(rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// determine the side of this ray with respect to the plane
		float t;
		const int side = (*rit).ComputeRayIntersection(axis, position, t);
	
		UpdateContributionsToPvs(*ray, true, side, pvsFront, pvsBack, pvsTotal);
		UpdateContributionsToPvs(*ray, false, side, pvsFront, pvsBack, pvsTotal);
	}

	//////////////
	//-- evaluate cost heuristics

	float pOverall = tData.mProbability;

	// we use spatial mid split => simplified computation
	pBack = pFront = pOverall * 0.5f;
	
	const float newCost = pvsBack * pBack + pvsFront * pFront;
	const float oldCost = (float)pvsTotal * pOverall + Limits::Small;
	
#ifdef GTPGTP_DEBUG
	Debug << "axis: " << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
	Debug << "p: " << pFront << " " << pBack << " " << pOverall << endl;
#endif

	return  (mCtDivCi + newCost) / oldCost;
}


void VspTree::UpdateContributionsToPvs(Intersectable *obj,
									   const int cf,
									   float &frontPvs,
									   float &backPvs,
									   float &totalPvs) const
{
	if (!obj) return;

	const float renderCost = mViewCellsManager->EvalRenderCost(obj);

	// object in no pvs => new
	if (!obj->Mailed() && !obj->Mailed(1) && !obj->Mailed(2))
	{
		totalPvs += renderCost;
	}

	// QUESTION matt: is it safe to assume that 
	// the object belongs to no pvs in this case?
	//if (cf == Ray::COINCIDENT) return;
	if (cf >= 0) // front pvs
	{
		if (!obj->Mailed() && !obj->Mailed(2))
		{
			frontPvs += renderCost;
		
			// already in back pvs => in both pvss
			if (obj->Mailed(1))
				obj->Mail(2);
			else
				obj->Mail();
		}
	}

	if (cf <= 0) // back pvs
	{
		if (!obj->Mailed(1) && !obj->Mailed(2))
		{
			backPvs += renderCost;
		
			// already in front pvs => in both pvss
			if (obj->Mailed())
				obj->Mail(2);
			else
				obj->Mail(1);
		}
	}
}


void VspTree::UpdateContributionsToPvs(BvhLeaf *leaf,
									   const int cf,
									   float &frontPvs,
									   float &backPvs,
									   float &totalPvs,
									   const bool countEntries) const
{
	if (!leaf) return;

	const float renderCost = countEntries ? 1 : 
		mHierarchyManager->mBvHierarchy->EvalAbsCost(leaf->mObjects);
	
	// leaf in no pvs => new
	if (!leaf->Mailed() && !leaf->Mailed(1) && !leaf->Mailed(2))
	{
		totalPvs += renderCost;
	}

	if (cf >= 0) // front pvs
	{
		if (!leaf->Mailed() && !leaf->Mailed(2))
		{
			frontPvs += renderCost;
	
			// already in back pvs => in both pvss
			if (leaf->Mailed(1))
				leaf->Mail(2);
			else
				leaf->Mail();
		}
	}

	if (cf <= 0) // back pvs
	{
		if (!leaf->Mailed(1) && !leaf->Mailed(2))
		{
			backPvs += renderCost;
		
			// already in front pvs => in both pvss
			if (leaf->Mailed())
				leaf->Mail(2);
			else
				leaf->Mail(1);
		}
	}
}


void VspTree::UpdateContributionsToPvs(KdLeaf *leaf,
									   const int cf,
									   float &frontPvs,
									   float &backPvs,
									   float &totalPvs) const
{
	if (!leaf) return;

	// the objects which are referenced in this and only this leaf
	const int contri = (int)(leaf->mObjects.size() - leaf->mMultipleObjects.size());
	
	// newly found leaf
	if (!leaf->Mailed() && !leaf->Mailed(1) && !leaf->Mailed(2))
	{
		totalPvs += contri;
	}

	// recursivly update contributions of yet unclassified objects
	ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();

	for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
	{	
		UpdateContributionsToPvs(*oit, cf, frontPvs, backPvs, totalPvs);
    }	
	
	if (cf >= 0) // front pvs
	{
		if (!leaf->Mailed() && !leaf->Mailed(2))
		{
			frontPvs += contri;
		
			// already in back pvs => in both pvss
			if (leaf->Mailed(1))
				leaf->Mail(2);
			else
				leaf->Mail();
		}
	}

	if (cf <= 0) // back pvs
	{
		if (!leaf->Mailed(1) && !leaf->Mailed(2))
		{
			backPvs += contri;
		
			// already in front pvs => in both pvss
			if (leaf->Mailed())
				leaf->Mail(2);
			else
				leaf->Mail(1);
		}
	}
}


void VspTree::CollectLeaves(vector<VspLeaf *> &leaves, 
							const bool onlyUnmailed,
							const int maxPvsSize) const
{
	stack<VspNode *> nodeStack;
	nodeStack.push(mRoot);

	while (!nodeStack.empty())
	{
		VspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
		{
			// test if this leaf is in valid view space
			VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);

			if (leaf->TreeValid() && 
				(!onlyUnmailed || !leaf->Mailed()) &&
				((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().EvalPvsCost() <= maxPvsSize)))
			{
				leaves.push_back(leaf);
			}
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);

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


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


VspNode *VspTree::GetRoot() const
{
	return mRoot;
}


void VspTree::EvaluateLeafStats(const VspTraversalData &data)
{
	// the node became a leaf -> evaluate stats for leafs
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);

	if (data.mPvs > mVspStats.maxPvs)
	{
		mVspStats.maxPvs = (int)data.mPvs;
	}

	mVspStats.pvs += (int)data.mPvs;

	if (data.mDepth < mVspStats.minDepth)
	{
		mVspStats.minDepth = data.mDepth;
	}
	
	if (data.mDepth >= mTermMaxDepth)
	{
        ++ mVspStats.maxDepthNodes;
	}

	// accumulate rays to compute rays /  leaf
	mVspStats.rayRefs += (int)data.mRays->size();

	if (data.mPvs < mTermMinPvs)
		++ mVspStats.minPvsNodes;

	if ((int)data.mRays->size() < mTermMinRays)
		++ mVspStats.minRaysNodes;

	if (data.GetAvgRayContribution() > mTermMaxRayContribution)
		++ mVspStats.maxRayContribNodes;

	if (data.mProbability <= mTermMinProbability)
		++ mVspStats.minProbabilityNodes;
	
	// accumulate depth to compute average depth
	mVspStats.accumDepth += data.mDepth;

	++ mCreatedViewCells;

#ifdef GTP_DEBUG
	Debug << "BSP stats: "
		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
		  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
		  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
		  << "#pvs: " << leaf->GetViewCell()->GetPvs().EvalPvsCost() << "), "
		  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
#endif
}


void VspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
{
	ViewCell::NewMail();
	CollectViewCells(mRoot, onlyValid, viewCells, true);
}


void VspTree::CollapseViewCells()
{
// TODO matt
#if HAS_TO_BE_REDONE
	stack<VspNode *> nodeStack;

	if (!mRoot)
		return;

	nodeStack.push(mRoot);
	
	while (!nodeStack.empty()) 
	{
		VspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
        {
			BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();

			if (!viewCell->GetValid())
			{
				BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
	
				ViewCellContainer leaves;
				mViewCellsTree->CollectLeaves(viewCell, leaves);

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

				for (it = leaves.begin(); it != it_end; ++ it)
				{
					VspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
					l->SetViewCell(GetOrCreateOutOfBoundsCell());
					++ mVspStats.invalidLeaves;
				}

				// add to unbounded view cell
				ViewCell *outOfBounds = GetOrCreateOutOfBoundsCell();
				outOfBounds->GetPvs().AddPvs(viewCell->GetPvs());
				DEL_PTR(viewCell);
			}
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
		
			nodeStack.push(interior->GetFront());
			nodeStack.push(interior->GetBack());
		}
	}

	Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
#endif
}


void VspTree::CollectRays(VssRayContainer &rays)
{
	vector<VspLeaf *> leaves;
	CollectLeaves(leaves);

	vector<VspLeaf *>::const_iterator lit, lit_end = leaves.end();

	for (lit = leaves.begin(); lit != lit_end; ++ lit)
	{
		VspLeaf *leaf = *lit;
		VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();

		for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
			rays.push_back(*rit);
	}
}


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


void VspTree::ValidateTree()
{
	if (!mRoot)
		return;

	mVspStats.invalidLeaves = 0;
	stack<VspNode *> nodeStack;

	nodeStack.push(mRoot);

	while (!nodeStack.empty()) 
	{
		VspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
		{
			VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);

			if (!leaf->GetViewCell()->GetValid())
				++ mVspStats.invalidLeaves;

			// validity flags don't match => repair
			if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
			{
				leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
				PropagateUpValidity(leaf);
			}
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
		
			nodeStack.push(interior->GetFront());
			nodeStack.push(interior->GetBack());
		}
	}

	Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
}



void VspTree::CollectViewCells(VspNode *root,
								  bool onlyValid,
								  ViewCellContainer &viewCells,
								  bool onlyUnmailed) const
{
	if (!root)
		return;

	stack<VspNode *> nodeStack;
	nodeStack.push(root);
	
	while (!nodeStack.empty()) 
	{
		VspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
		{
			if (!onlyValid || node->TreeValid())
			{
				ViewCellLeaf *leafVc = dynamic_cast<VspLeaf *>(node)->GetViewCell();

				ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
						
				if (!onlyUnmailed || !viewCell->Mailed()) 
				{
					viewCell->Mail();
					viewCells.push_back(viewCell);
				}
			}
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
		
			nodeStack.push(interior->GetFront());
			nodeStack.push(interior->GetBack());
		}
	}
}


int VspTree::FindNeighbors(VspLeaf *n,
						   vector<VspLeaf *> &neighbors,
						   const bool onlyUnmailed) const
{
	stack<VspNode *> nodeStack;
	nodeStack.push(mRoot);

	const AxisAlignedBox3 box = GetBoundingBox(n);

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

		if (node->IsLeaf())
		{
			VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);

			if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
				neighbors.push_back(leaf);
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
			
			if (interior->GetPosition() > box.Max(interior->GetAxis()))
				nodeStack.push(interior->GetBack());
			else
			{
				if (interior->GetPosition() < box.Min(interior->GetAxis()))
					nodeStack.push(interior->GetFront());
				else
				{
					// random decision
					nodeStack.push(interior->GetBack());
					nodeStack.push(interior->GetFront());
				}
			}
		}
	}

	return (int)neighbors.size();
}


// Find random neighbor which was not mailed
VspLeaf *VspTree::GetRandomLeaf(const Plane3 &plane)
{
	stack<VspNode *> nodeStack;
	nodeStack.push(mRoot);
  
	int mask = rand();
  
	while (!nodeStack.empty()) 
	{
		VspNode *node = nodeStack.top();
		
		nodeStack.pop();
		
		if (node->IsLeaf()) 
		{
			return dynamic_cast<VspLeaf *>(node);
		} 
		else 
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
			VspNode *next;
			
			if (GetBoundingBox(interior->GetBack()).Side(plane) < 0)
			{
				next = interior->GetFront();
			}
            else
			{
				if (GetBoundingBox(interior->GetFront()).Side(plane) < 0)
				{
					next = interior->GetBack();
				}
				else 
				{
					// random decision
					if (mask & 1)
						next = interior->GetBack();
					else
						next = interior->GetFront();
					mask = mask >> 1;
				}
			}
			
			nodeStack.push(next);
		}
	}
  
	return NULL;
}


VspLeaf *VspTree::GetRandomLeaf(const bool onlyUnmailed)
{
	stack<VspNode *> nodeStack;

	nodeStack.push(mRoot);

	int mask = rand();

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

		if (node->IsLeaf())
		{
			if ( (!onlyUnmailed || !node->Mailed()) )
				return dynamic_cast<VspLeaf *>(node);
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);

			// random decision
			if (mask & 1)
				nodeStack.push(interior->GetBack());
			else
				nodeStack.push(interior->GetFront());

			mask = mask >> 1;
		}
	}

	return NULL;
}


float VspTree::EvalPvsCost(const RayInfoContainer &rays) const
{
	float pvsCost = 0;
	
	Intersectable::NewMail();
	KdNode::NewMail();

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

	for (rit = rays.begin(); rit != rays.end(); ++ rit)
	{
		VssRay *ray = (*rit).mRay;
		
		pvsCost += EvalContributionToPvs(*ray, true);

#if COUNT_ORIGIN_OBJECTS
		pvsCost += EvalContributionToPvs(*ray, false);
#endif
	}
	
	return pvsCost;
}


int VspTree::EvalPvsEntriesContribution(const VssRay &ray, 
										const bool isTermination) const

{
		Intersectable *obj; 
		Vector3 pt;
		KdNode *node;

		ray.GetSampleData(isTermination, pt, &obj, &node);
		if (!obj) return 0;

		switch(mHierarchyManager->GetObjectSpaceSubdivisionType())
		{
		case HierarchyManager::NO_OBJ_SUBDIV:
		{
			if (!obj->Mailed())
			{
				obj->Mail();
				return 1;
			}
			
			return 0;
		}

		case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);
			if (!leaf->Mailed())
			{
				leaf->Mail();
				return 1;
			}
			
			return 0;
		}
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *bvhleaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);

			if (!bvhleaf->Mailed())
			{
				bvhleaf->Mail();
				return 1;
			}
			
			return 0;
		}
	default:
		break;
	}

	return 0;
}


int VspTree::EvalPvsEntriesSize(const RayInfoContainer &rays) const
{
	int pvsSize = 0;

	Intersectable::NewMail();
	KdNode::NewMail();

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

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

		pvsSize += EvalPvsEntriesContribution(*ray, true);

#if COUNT_ORIGIN_OBJECTS
		pvsSize += EvalPvsEntriesContribution(*ray, false);
#endif
	}

	return pvsSize;
}


float VspTree::GetEpsilon() const
{
	return mEpsilon;
}


int VspTree::CastLineSegment(const Vector3 &origin,
							 const Vector3 &termination,
                             ViewCellContainer &viewcells,
							 const bool useMailboxing)
{
	int hits = 0;

	float mint = 0.0f, maxt = 1.0f;
	const Vector3 dir = termination - origin;

	stack<LineTraversalData> tStack;

	Vector3 entp = origin;
	Vector3 extp = termination;

	VspNode *node = mRoot;
	VspNode *farChild;

	float position;
	int axis;

	while (1)
	{
		if (!node->IsLeaf())
		{
			VspInterior *in = dynamic_cast<VspInterior *>(node);
			position = in->GetPosition();
			axis = in->GetAxis();

			if (entp[axis] <= position)
			{
				if (extp[axis] <= position)
				{
					node = in->GetBack();
					// cases N1,N2,N3,P5,Z2,Z3
					continue;
				} else
				{
					// case N4
					node = in->GetBack();
					farChild = in->GetFront();
				}
			}
			else
			{
				if (position <= extp[axis])
				{
					node = in->GetFront();
					// cases P1,P2,P3,N5,Z1
					continue;
				}
				else
				{
					node = in->GetFront();
					farChild = in->GetBack();
					// case P4
				}
			}

 			// $$ modification 3.5.2004 - hints from Kamil Ghais
			// case N4 or P4
			const float tdist = (position - origin[axis]) / dir[axis];
			tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO

			extp = origin + dir * tdist;
			maxt = tdist;
		}
		else
		{
			// compute intersection with all objects in this leaf
			VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
			ViewCell *viewCell;
			if (0)
				viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
			else
				viewCell = leaf->GetViewCell();

			// don't have to mail if each view cell belongs to exactly one leaf
			if (!useMailboxing || !viewCell->Mailed())
			{
				if (useMailboxing)
					viewCell->Mail();

				viewcells.push_back(viewCell);
				++ hits;
			}

			// get the next node from the stack
			if (tStack.empty())
				break;

			entp = extp;
			mint = maxt;
			
			LineTraversalData &s  = tStack.top();
			node = s.mNode;
			extp = s.mExitPoint;
			maxt = s.mMaxT;

			tStack.pop();
		}
	}

	return hits;
}


int VspTree::CastRay(Ray &ray)
{
	int hits = 0;

	stack<LineTraversalData> tStack;
	const Vector3 dir = ray.GetDir();

	float maxt, mint;

	if (!mBoundingBox.GetRaySegment(ray, mint, maxt))
		return 0;

	Intersectable::NewMail();
	ViewCell::NewMail();

	Vector3 entp = ray.Extrap(mint);
	Vector3 extp = ray.Extrap(maxt);

	const Vector3 origin = entp;

	VspNode *node = mRoot;
	VspNode *farChild = NULL;

	float position;
	int axis;

	while (1)
	{
		if (!node->IsLeaf())
		{
			VspInterior *in = dynamic_cast<VspInterior *>(node);
			position = in->GetPosition();
			axis = in->GetAxis();

			if (entp[axis] <= position)
			{
				if (extp[axis] <= position)
				{
					node = in->GetBack();
					// cases N1,N2,N3,P5,Z2,Z3
					continue;
				} 
				else
				{
					// case N4
					node = in->GetBack();
					farChild = in->GetFront();
				}
			}
			else
			{
				if (position <= extp[axis])
				{
					node = in->GetFront();
					// cases P1,P2,P3,N5,Z1
					continue;
				}
				else
				{
					node = in->GetFront();
					farChild = in->GetBack();
					// case P4
				}
			}

 			// $$ modification 3.5.2004 - hints from Kamil Ghais
			// case N4 or P4
			const float tdist = (position - origin[axis]) / dir[axis];
			tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
			extp = origin + dir * tdist;
			maxt = tdist;
		}
		else
		{
			// compute intersection with all objects in this leaf
			VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
			ViewCell *vc = leaf->GetViewCell();

			if (!vc->Mailed())
			{
				vc->Mail();
				// todo: add view cells to ray
				++ hits;
			}

			// get the next node from the stack
			if (tStack.empty())
				break;

			entp = extp;
			mint = maxt;
			
			LineTraversalData &s  = tStack.top();
			node = s.mNode;
			extp = s.mExitPoint;
			maxt = s.mMaxT;
			tStack.pop();
		}
	}

	return hits;
}


ViewCell *VspTree::GetViewCell(const Vector3 &point, const bool active)
{
	if (mRoot == NULL)
		return NULL;

	stack<VspNode *> nodeStack;
	nodeStack.push(mRoot);
  
	ViewCellLeaf *viewcell = NULL;
  
	while (!nodeStack.empty())  
	{
		VspNode *node = nodeStack.top();
		nodeStack.pop();
	
		if (node->IsLeaf()) 
		{
			/*const AxisAlignedBox3 box = GetBoundingBox(dynamic_cast<VspLeaf *>(node));
			if (!box.IsInside(point))
				cerr << "error, point " << point << " should be in view cell " << box << endl;
			*/	
			viewcell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
			break;
		} 
		else 	
		{	
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
     
			// random decision
			if (interior->GetPosition() - point[interior->GetAxis()] < 0)
			{
				nodeStack.push(interior->GetFront());
			}
			else
			{
				nodeStack.push(interior->GetBack());
			}
		}
	}
  
	if (active)
	{
		return mViewCellsTree->GetActiveViewCell(viewcell);
	}
	else
	{
		return viewcell;
	}
}


bool VspTree::ViewPointValid(const Vector3 &viewPoint) const
{
	VspNode *node = mRoot;

	while (1)
	{
		// early exit
		if (node->TreeValid())
			return true;

		if (node->IsLeaf())
			return false;
			
		VspInterior *in = dynamic_cast<VspInterior *>(node);
					
		if (in->GetPosition() - viewPoint[in->GetAxis()] <= 0) 
		{
			node = in->GetBack();
		} 
		else
		{
			node = in->GetFront();
		}
	}

	// should never come here
	return false;
}


void VspTree::PropagateUpValidity(VspNode *node)
{
	const bool isValid = node->TreeValid();

	// propagative up invalid flag until only invalid nodes exist over this node
	if (!isValid)
	{
		while (!node->IsRoot() && node->GetParent()->TreeValid())
		{
			node = node->GetParent();
			node->SetTreeValid(false);
		}
	}
	else
	{
		// propagative up valid flag until one of the subtrees is invalid
		while (!node->IsRoot() && !node->TreeValid())
		{
            node = node->GetParent();
			VspInterior *interior = dynamic_cast<VspInterior *>(node);
			
			// the parent is valid iff both leaves are valid
			node->SetTreeValid(interior->GetBack()->TreeValid() && 
							   interior->GetFront()->TreeValid());
		}
	}
}


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

	return true;
}


void VspTree::ExportNode(VspNode *node, OUT_STREAM &stream)
{
	if (node->IsLeaf())
	{
		VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
		ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());

		int id = -1;
		if (viewCell != mOutOfBoundsCell)
			id = viewCell->GetId();

		stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
	}
	else
	{	
		VspInterior *interior = dynamic_cast<VspInterior *>(node);
	
		AxisAlignedPlane plane = interior->GetPlane();
		stream << "<Interior plane=\"" << plane.mPosition << " " 
			   << plane.mAxis << "\">" << endl;

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

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


int VspTree::SplitRays(const AxisAlignedPlane &plane,
					   RayInfoContainer &rays,
					   RayInfoContainer &frontRays,
					   RayInfoContainer &backRays) const
{
	int splits = 0;

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

	for (rit = rays.begin(); rit != rit_end; ++ rit)
	{
		RayInfo bRay = *rit;
		
		VssRay *ray = bRay.mRay;
		float t;

		// get classification and receive new t
		// test if start point behind or in front of plane
		const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
			
		if (side == 0)
		{
			++ splits;

			if (ray->HasPosDir(plane.mAxis))
			{
				backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
				frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
			}
			else
			{
				frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
				backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
			}
		}
		else if (side == 1)
		{
			frontRays.push_back(bRay);
		}
		else
		{
			backRays.push_back(bRay);
		}
	}

	return splits;
}


AxisAlignedBox3 VspTree::GetBoundingBox(VspNode *node) const
{
	if (!node->GetParent())
		return mBoundingBox;

	if (!node->IsLeaf())
	{
		return (dynamic_cast<VspInterior *>(node))->GetBoundingBox();		
	}

	VspInterior *parent = dynamic_cast<VspInterior *>(node->GetParent());

	AxisAlignedBox3 box(parent->GetBoundingBox());

	if (parent->GetFront() == node)
		box.SetMin(parent->GetAxis(), parent->GetPosition());
    else
		box.SetMax(parent->GetAxis(), parent->GetPosition());

	return box;
}


int VspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box,
									 ViewCellContainer &viewCells) const
{
	stack<VspNode *> nodeStack;
 
	ViewCell::NewMail();

	nodeStack.push(mRoot);
	
	while (!nodeStack.empty()) 
	{
		VspNode *node = nodeStack.top();
		nodeStack.pop();

		const AxisAlignedBox3 bbox = GetBoundingBox(node);

		if (Overlap(bbox, box)) {
		  if (node->IsLeaf())
			{
			  VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
			  
			  if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
				{
				  leaf->GetViewCell()->Mail();
				  viewCells.push_back(leaf->GetViewCell());
				}
			}
		  else 
			{
			  VspInterior *interior = dynamic_cast<VspInterior *>(node);
			  
			  VspNode *first = interior->GetFront();
			  VspNode *second = interior->GetBack();
			  
			  nodeStack.push(first);
			  nodeStack.push(second);
			}
		}
		// default: cull
	}

	return (int)viewCells.size();
}


void VspTree::PreprocessRays(const VssRayContainer &sampleRays,
							 RayInfoContainer &rays)
{
	VssRayContainer::const_iterator rit, rit_end = sampleRays.end();

	long startTime = GetTime();

	cout << "storing rays ... ";

	Intersectable::NewMail();

	//-- store rays and objects
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = *rit;
		float minT, maxT;
		static Ray hray;

		hray.Init(*ray);
		
		// TODO: not very efficient to implictly cast between rays types
		if (GetBoundingBox().GetRaySegment(hray, minT, maxT))
		{
			float len = ray->Length();

			if (!len)
				len = Limits::Small;

			rays.push_back(RayInfo(ray, minT / len, maxT / len));
		}
	}

	cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
}


void VspTree::GetViewCells(const VssRay &ray, ViewCellContainer &viewCells)
{
	static Ray hray;
	hray.Init(ray);
	
	float tmin = 0, tmax = 1.0;

	if (!mBoundingBox.GetRaySegment(hray, tmin, tmax) || (tmin > tmax))
		return;

	const Vector3 origin = hray.Extrap(tmin);
	const Vector3 termination = hray.Extrap(tmax);

	// view cells were not precomputed
	// don't mail because we need mailboxing for something else
	CastLineSegment(origin, termination, viewCells, false);
}


void VspTree::Initialise(const VssRayContainer &rays,
						 AxisAlignedBox3 *forcedBoundingBox)
{
	ComputeBoundingBox(rays, forcedBoundingBox);

	VspLeaf *leaf = new VspLeaf();
	mRoot = leaf;

	VspViewCell *viewCell = new VspViewCell();
    leaf->SetViewCell(viewCell);

	// set view cell values
	viewCell->mLeaves.push_back(leaf);
	viewCell->SetVolume(mBoundingBox.GetVolume());

	leaf->mProbability = mBoundingBox.GetVolume();
}


void VspTree::PrepareConstruction(SplitQueue &tQueue,
								  const VssRayContainer &sampleRays,
								  RayInfoContainer &rays)
{	
	mVspStats.Reset();
	mVspStats.Start();
	mVspStats.nodes = 1;

	// store pointer to this tree
	VspSubdivisionCandidate::sVspTree = this;
	
	// initialise termination criteria
	mTermMinProbability *= mBoundingBox.GetVolume();
	
	// get clipped rays
	PreprocessRays(sampleRays, rays);

	/// collect pvs from rays
	const float pvsCost = EvalPvsCost(rays);
	
	// root and bounding box were already constructed
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(mRoot);

	//////////
	//-- prepare view space partition

	const float prop = mBoundingBox.GetVolume();
	
	// first vsp traversal data
	VspTraversalData vData(leaf, 0, &rays, pvsCost, prop, mBoundingBox);

#if WORK_WITH_VIEWCELL_PVS
	// add first view cell to all the objects view cell pvs
	ObjectPvsEntries::const_iterator oit, 
		oit_end = leaf->GetViewCell()->GetPvs().mEntries.end();

	for (oit = leaf->GetViewCell()->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = (*oit).first;
		obj->mViewCellPvs.AddSample(leaf->GetViewCell(), 1);
	}
#endif

	mTotalCost = vData.mCorrectedRenderCost = vData.mRenderCost = pvsCost;
	mPvsEntries = vData.mCorrectedPvs = vData.mPvs = EvalPvsEntriesSize(rays);

	//////////////
	//-- create the first split candidate

	VspSubdivisionCandidate *splitCandidate = new VspSubdivisionCandidate(vData);
    EvalSubdivisionCandidate(*splitCandidate);
	leaf->SetSubdivisionCandidate(splitCandidate);

	EvalSubdivisionStats(*splitCandidate);

	tQueue.Push(splitCandidate);
}


void VspTree::CollectDirtyCandidate(const VssRay &ray, 
									const bool isTermination,
									vector<SubdivisionCandidate *> &dirtyList,
									const bool onlyUnmailed) const
{

	Intersectable *obj;
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);
	
	if (!obj) return;
	
	SubdivisionCandidate *candidate = NULL;
		
	switch (mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
	case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);

			if (!leaf->Mailed())
			{
				leaf->Mail();
				candidate = leaf->mSubdivisionCandidate;
			}
			break;
		}
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *leaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);

			if (!leaf->Mailed())
			{
				leaf->Mail();
				candidate = leaf->GetSubdivisionCandidate();
			}
			break;
		}
	default:
		cerr << "not implemented yet" << endl;
		candidate = NULL;
		break;
	}

	// is this leaf still a split candidate?
	if (candidate && (!onlyUnmailed || !candidate->Mailed()))
	{
		candidate->Mail();
		candidate->SetDirty(true);
		dirtyList.push_back(candidate);
	}
}


void VspTree::CollectDirtyCandidates(VspSubdivisionCandidate *sc, 
									 vector<SubdivisionCandidate *> &dirtyList,
									 const bool onlyUnmailed)
{
	VspTraversalData &tData = sc->mParentData;
	VspLeaf *node = tData.mNode;
	
	KdLeaf::NewMail();
	Intersectable::NewMail();
	
	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();

	// add all kd nodes seen by the rays
	for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
		
		CollectDirtyCandidate(*ray, true, dirtyList, onlyUnmailed);
        CollectDirtyCandidate(*ray, false, dirtyList, onlyUnmailed);
	}
}


int VspTree::EvalMaxEventContribution(const VssRay &ray, 
									  const bool isTermination) const
{
	Intersectable *obj;
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);

	if (!obj) 
		return 0;

	int pvs = 0;

	switch (mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
	case HierarchyManager::NO_OBJ_SUBDIV:
		{
			if (-- obj->mCounter == 0)
				++ pvs;
			break;
		}
	case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);

			// add contributions of the kd nodes
			pvs += EvalMaxEventContribution(leaf);
			break;
		}
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *leaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);

			if (-- leaf->mCounter == 0)
				pvs += (int)leaf->mObjects.size();
			break;
		}
	default:
		break;
	}

	return pvs;
}


float VspTree::PrepareHeuristics(const VssRay &ray, const bool isTermination)
{
	float pvsSize = 0;
	
	Intersectable *obj;
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);

	if (!obj)
		return 0;

	switch (mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
	case HierarchyManager::NO_OBJ_SUBDIV:
		{
			if (!obj->Mailed())
			{
				obj->Mail();
				obj->mCounter = 0;
				++ pvsSize;
			}

			++ obj->mCounter;	
			break;
		}
	case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);
			pvsSize += PrepareHeuristics(leaf);	
			break;
		}
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *leaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);

			if (!leaf->Mailed())
			{
				leaf->Mail();
				leaf->mCounter = 0;
				pvsSize += (int)leaf->mObjects.size();
			}

			++ leaf->mCounter;	
			break;
		}
	default:
		break;
	}

	return pvsSize;
}


int VspTree::EvalMinEventContribution(const VssRay &ray, 
									  const bool isTermination) const
{
	Intersectable *obj;
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);

	if (!obj) return 0;

	int pvs = 0;

	switch (mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
	case HierarchyManager::NO_OBJ_SUBDIV:
		{
			if (!obj->Mailed())
			{
				obj->Mail();
				++ pvs;
			}
			break;
		}
	case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);
			// add contributions of the kd nodes
			pvs += EvalMinEventContribution(leaf);				
			break;
		}
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *leaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);
			if (!leaf->Mailed())
			{
				leaf->Mail();
				pvs += (int)leaf->mObjects.size();
			}
			break;
		}
	default:
		break;
	}

	return pvs;
}


void VspTree::UpdateContributionsToPvs(const VssRay &ray,
									   const bool isTermination,
									   const int cf,
									   float &frontPvs,
									   float &backPvs,
									   float &totalPvs) const
{
	Intersectable *obj;
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);

	if (!obj) return;

	switch (mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
		case HierarchyManager::NO_OBJ_SUBDIV:
		{
			// find front and back pvs for origing and termination object
			UpdateContributionsToPvs(obj, cf, frontPvs, backPvs, totalPvs);
			break;
		}
		case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);
			UpdateContributionsToPvs(leaf, cf, frontPvs, backPvs, totalPvs);
			break;
		}
		case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *leaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);
			UpdateContributionsToPvs(leaf, cf, frontPvs, backPvs, totalPvs);
			break;
		}
	}
}


void VspTree::UpdatePvsEntriesContribution(const VssRay &ray,
										   const bool isTermination,
										   const int cf,
										   float &pvsFront,
										   float &pvsBack,
										   float &totalPvs) const
{
	Intersectable *obj;
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);
	if (!obj) return;

	switch (mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
	case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		// TODO
		break;
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *leaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);
			UpdateContributionsToPvs(leaf, cf, pvsFront, pvsBack, totalPvs, true);
			break;
		}
	default:
		UpdateContributionsToPvs(obj, cf, pvsFront, pvsBack, totalPvs);	
		break;
	}
}


int VspTree::EvalContributionToPvs(const VssRay &ray, const bool isTermination) const
{	
	Intersectable *obj; 
	Vector3 pt;
	KdNode *node;

	ray.GetSampleData(isTermination, pt, &obj, &node);

	if (!obj) return 0;

	int pvs = 0;

	switch(mHierarchyManager->GetObjectSpaceSubdivisionType())
	{
	case HierarchyManager::NO_OBJ_SUBDIV:
		{
			if (!obj->Mailed())
			{
				obj->Mail();
				++ pvs;
			}
			break;
		}
	case HierarchyManager::KD_BASED_OBJ_SUBDIV:
		{
			KdLeaf *leaf = mHierarchyManager->mOspTree->GetLeaf(pt, node);
			pvs += EvalContributionToPvs(leaf);
			break;
		}
	case HierarchyManager::BV_BASED_OBJ_SUBDIV:
		{
			BvhLeaf *bvhleaf = mHierarchyManager->mBvHierarchy->GetLeaf(obj);

			if (!bvhleaf->Mailed())
			{
				bvhleaf->Mail();
				pvs += (int)bvhleaf->mObjects.size();
			}
			break;
		}
	default:
		break;
	}

	return pvs;
}


int VspTree::EvalContributionToPvs(KdLeaf *leaf) const
{
	if (leaf->Mailed()) // leaf already mailed
		return 0;
	
	leaf->Mail();

	// this is the pvs which is uniquely part of this kd leaf
	int pvs = (int)(leaf->mObjects.size() - leaf->mMultipleObjects.size());

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

	for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		if (!obj->Mailed())
		{
			obj->Mail();
			++ pvs;
		}
	}

	return pvs;
}


VspNode *VspTree::SubdivideAndCopy(SplitQueue &tQueue,
								   SubdivisionCandidate *splitCandidate)
{
	// todo remove dynamic cast
	VspSubdivisionCandidate *sc = dynamic_cast<VspSubdivisionCandidate *>(splitCandidate);

	VspTraversalData &tData = sc->mParentData;
	VspNode *newNode = tData.mNode;
	VspNode *oldNode = (VspNode *)splitCandidate->mEvaluationHack;

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

		VspTraversalData tFrontData;
		VspTraversalData tBackData;
		
		VspInterior *oldInterior = dynamic_cast<VspInterior *>(oldNode);

		// create new interior node and two leaf node
		const AxisAlignedPlane splitPlane = oldInterior->GetPlane();
		sc->mSplitPlane = splitPlane;
	
		// evaluate the changes in render cost and pvs entries
		EvalSubdivisionCandidate(*sc, false);

		newNode = SubdivideNode(*sc, tFrontData, tBackData);
	
		//oldNode->mRenderCostDecr += sc->GetRenderCostDecrease();
		//oldNode->mPvsEntriesIncr += sc->GetPvsEntriesIncr();

		//oldNode->mRenderCostDecr = sc->GetRenderCostDecrease();
		//oldNode->mPvsEntriesIncr = sc->GetPvsEntriesIncr();

		/////////////
		//-- evaluate new split candidates for global greedy cost heuristics

		VspSubdivisionCandidate *frontCandidate = new VspSubdivisionCandidate(tFrontData);
		VspSubdivisionCandidate *backCandidate = new VspSubdivisionCandidate(tBackData);

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

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

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

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

		// note: leaf is not destroyed because it is needed to collect 
		// dirty candidates in hierarchy manager
	}

	if (newNode->IsLeaf()) // subdivision terminated
	{
		// detach subdivision candidate: this leaf is no candidate for splitting anymore
		tData.mNode->SetSubdivisionCandidate(NULL); 
		// detach node so it won't get deleted
		tData.mNode = NULL;
	}

	return newNode;
}


int VspTree::CompressObjects(VspLeaf *leaf)
{
	bool compressed = true;

	while (compressed)
	{
		BvhNode::NewMail(2);

		ObjectPvsIterator oit = leaf->GetViewCell()->GetPvs().GetIterator();
		vector<BvhNode *> parents;
		ObjectPvs newPvs;

		compressed = false;

		while (oit.HasMoreEntries())
		{
			const ObjectPvsEntry &entry = oit.Next();
			BvhNode *obj = dynamic_cast<BvhNode *>(entry.mObject);

			if (!obj->IsRoot())
			{
				BvhNode *parent = obj->GetParent();

				if (!parent->Mailed())
				{
					if (parent->Mailed(1))
						cout << "error!!" << endl;
					parent->Mail();
				}
				else
				{
					// sibling was already found => 
					// this entry can be exchanged by the parent
					parents.push_back(parent);
					parent->Mail(1);
		
					compressed = true;
				}
			}
		}
		
		oit = leaf->GetViewCell()->GetPvs().GetIterator();

		while (oit.HasMoreEntries())
		{
			const ObjectPvsEntry &entry = oit.Next();
			BvhNode *obj = dynamic_cast<BvhNode *>(entry.mObject);

			if (!obj->IsRoot())
			{
				BvhNode *parent = obj->GetParent();

				// add only entries that cannot be exchaned with the parent
				if (!parent->Mailed(1))
				{
					newPvs.AddSampleDirty(obj, entry.mData.mSumPdf);
				}
			}
		}

		// add parents
		vector<BvhNode *>::const_iterator bit, bit_end = parents.end();

		for (bit = parents.begin(); bit != bit_end; ++ bit)
		{
			newPvs.AddSampleDirty(*bit, 1);
		}

		//cout << "size " << newPvs.GetSize() << endl;
		leaf->GetViewCell()->SetPvs(newPvs);
	}

	return leaf->GetViewCell()->GetPvs().GetSize();
}


int VspTree::CompressObjects()
{
	vector<VspLeaf *> leaves;
	CollectLeaves(leaves);

	int numEntries = 0;

	vector<VspLeaf *>::const_iterator lit, lit_end = leaves.end();

	for (lit = leaves.begin(); lit != lit_end; ++ lit)
	{
		numEntries += CompressObjects(*lit);
	}

	return numEntries;
}

}