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

#include "ViewCell.h"
#include "Plane3.h"
#include "VspOspTree.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 "KdIntersectable.h"


namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0


//-- static members

VspTree *VspTree::VspSplitCandidate::sVspTree = NULL;
OspTree *OspTree::OspSplitCandidate::sOspTree = NULL;

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

int VspNode::sMailId = 1;

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

	return (float)pvs;
}


static bool ViewCellHasMultipleReferences(Intersectable *obj, 
										  ViewCell *vc, 
										  bool checkOnlyMailed)
{
	MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);
//return false;
	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 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 << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits[0] + splits[1] + splits[2] << endl;

	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 << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
	//app << "#N_PVS: " << pvs << endl;

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

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


void OspTreeStatistics::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[0] + splits[1] + splits[2] << endl;

	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_PMINOBJECTREFLEAVES  ( Percentage of leaves with minimal number of object ref)\n" 
		//<< minObjectNodes * 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 << "#N_PVS: " << pvs << endl;

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


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


VspNode::VspNode(): 
mParent(NULL), mTreeValid(true), mTimeStamp(0)
{}


VspNode::VspNode(VspInterior *parent): 
mParent(parent), mTreeValid(true)
{}


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), mPvs(NULL)
{
}


VspLeaf::~VspLeaf()
{
	DEL_PTR(mPvs);
	CLEAR_CONTAINER(mVssRays);
}


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


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


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



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

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),
mOspTree(NULL)
{
	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("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);

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

	//-- partition criteria
	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.epsilon", mEpsilon);
	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()->GetIntValue("VspTree.maxTests", mMaxTests);
	Environment::GetSingleton()->GetFloatValue("VspTree.maxStaticMemory", mMaxMemory);

	Environment::GetSingleton()->GetBoolValue("VspTree.useCostHeuristics", mUseCostHeuristics);
	Environment::GetSingleton()->GetBoolValue("VspTree.simulateOctree", mCirculatingAxis);
	
	Environment::GetSingleton()->GetBoolValue("VspTree.useKdPvsForHeuristics", mUseKdPvsForHeuristics);

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

	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.minBand", mMinBand);
	Environment::GetSingleton()->GetFloatValue("VspTree.Construction.maxBand", mMaxBand);
	

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

	if (!mUseKdPvsForHeuristics)
		Debug << "pvs heuristics: per object" << endl;
	else
		Debug << "pvs heuristics: per kd node" << endl;

	mLocalSplitCandidates = new vector<SortableEntry>;

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


void VspTree::ComputeBoundingBox(const VssRayContainer &rays,
								 AxisAlignedBox3 *forcedBoundingBox) 
{	
	if (forcedBoundingBox)
	{
		mBoundingBox = *forcedBoundingBox;
	}
	else // compute vsp tree bounding box
	{
		mBoundingBox.Initialize();
		VssRayContainer::const_iterator rit, rit_end = rays.end();

		//-- compute bounding box
        for (rit = rays.begin(); rit != rit_end; ++ rit)
		{
			VssRay *ray = *rit;

			// compute bounding box of view space
			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;
}


// TODO: return memory usage in MB
float VspTree::GetMemUsage() const
{
	return (float)
		 (sizeof(VspTree) + 
		  mVspStats.Leaves() * sizeof(VspLeaf) + 
		  mCreatedViewCells * sizeof(VspViewCell) +
		  mVspStats.pvs * sizeof(PvsData) +
		  mVspStats.Interior() * sizeof(VspInterior) +
		  mVspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
}


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

	if (0 && localTerminationCriteriaMet)
	{
		Debug << "********local termination *********" << 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;
	}

	return localTerminationCriteriaMet;
		
}


inline bool VspTree::GlobalTerminationCriteriaMet(const VspTraversalData &data) const
{
	const bool terminationCriteriaMet = (
		// mOutOfMemory || 
		(mVspStats.Leaves() >= mMaxViewCells) || 
        (mGlobalCostMisses >= mTermGlobalCostMissTolerance) 
		);

	if (0 && terminationCriteriaMet)
	{
		Debug << "********* terminationCriteriaMet *********" << endl;
		Debug << "cost misses: " << mGlobalCostMisses << " " << mTermGlobalCostMissTolerance << endl;
		Debug << "leaves: " << mVspStats.Leaves() << " " <<  mMaxViewCells << endl;
	}
	return terminationCriteriaMet;
}


void VspTree::CreateViewCell(VspTraversalData &tData, const bool updatePvs)
{
	//-- create new view cell
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(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.CountObjectsInPvs(), pvs.GetSize());

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


	//-- store additional info
	if (mStoreRays)
	{
		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);
		}
	}
		
	// set view cell values
	viewCell->mLeaf = leaf;

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


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


VspNode *VspTree::Subdivide(SplitQueue &tQueue,
							SplitCandidate *splitCandidate,
							const bool globalCriteriaMet)
{
	// todo remove dynamic cast
	VspSplitCandidate *sc = dynamic_cast<VspSplitCandidate *>(splitCandidate);
	VspTraversalData &tData = sc->mParentData;

	VspNode *newNode = tData.mNode;


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

		//-- continue subdivision
		
		// create new interior node and two leaf node
		const AxisAlignedPlane splitPlane = sc->mSplitPlane;
		newNode = SubdivideNode(splitPlane, 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();
		mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;

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

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

		EvalSplitCandidate(*frontCandidate);
		EvalSplitCandidate(*backCandidate);

		tQueue.Push(frontCandidate);
		tQueue.Push(backCandidate);
		
		// delete old view cell
		delete tData.mNode->mViewCell;

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

	if (newNode->IsLeaf()) // subdivision terminated
	{
		// view cell is created during subdivision
		//CreateViewCell(tData); 

		VspLeaf *leaf = dynamic_cast<VspLeaf *>(newNode);
		ViewCell *viewCell = leaf->GetViewCell();

		int conSamp = 0;
		float sampCon = 0.0f;

#if 0
		//-- store pvs optained from rays

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

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

		mVspStats.contributingSamples += conSamp;
		mVspStats.sampleContributions += (int)sampCon;
#endif
		//-- store additional info
		if (mStoreRays)
		{
			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);
			}
		}

		// finally evaluate statistics for this leaf
		EvaluateLeafStats(tData);
	}

	//-- cleanup
	tData.Clear();
	
	return newNode;
}


void VspTree::EvalSplitCandidate(VspSplitCandidate &splitCandidate)
{
	float frontProb;
	float backProb;
	
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(splitCandidate.mParentData.mNode);
	
	// compute locally best split plane
	const bool success = 
		SelectSplitPlane(splitCandidate.mParentData, splitCandidate.mSplitPlane, frontProb, backProb);

	float oldRenderCost;

	// compute global decrease in render cost
	const float renderCostDecr = EvalRenderCostDecrease(splitCandidate.mSplitPlane, 
														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
	
	splitCandidate.SetPriority(priority);

	// max cost threshold violated?
	splitCandidate.mMaxCostMisses = 
		success ? splitCandidate.mParentData.mMaxCostMisses : splitCandidate.mParentData.mMaxCostMisses + 1;
	//Debug << "p: " << tData.mNode << " depth: " << tData.mDepth << endl;
}


VspInterior *VspTree::SubdivideNode(const AxisAlignedPlane &splitPlane,
									VspTraversalData &tData,
									VspTraversalData &frontData,
									VspTraversalData &backData)
{
	VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
	
	//-- the front and back traversal data is filled with the new 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);
	
	//Debug << "f: " << frontData.mRays->size() << " b: " << backData.mRays->size() << "d: " << tData.mRays->size() << endl;

	//-- compute pvs
	frontData.mPvs = EvalPvsSize(*frontData.mRays);
	backData.mPvs = EvalPvsSize(*backData.mRays);

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

	
    ///////////////////////////////////////////
	// subdivide further

	// store maximal and minimal depth
	if (tData.mDepth > mVspStats.maxDepth)
	{
		Debug << "max depth increases to " << tData.mDepth << " at " << mVspStats.Leaves() << " leaves" << endl;
		mVspStats.maxDepth = tData.mDepth;
	}

	// two more leaves
	mVspStats.nodes += 2;
    
	VspInterior *interior = new VspInterior(splitPlane);

#ifdef _DEBUG
	Debug << interior << endl;
#endif+


	//-- create front and back leaf

	VspInterior *parent = leaf->GetParent();

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

		// remove "parent" view cell from pvs of all objects (traverse trough rays)
		RemoveParentViewCellReferences(tData.mNode->GetViewCell());
	}
	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;

	// explicitely create front and back view cell
	CreateViewCell(frontData, false);
	CreateViewCell(backData, false);


#if WORK_WITH_VIEWCELL_PVS
	// create front and back view cell
	// add front and back view cell to "Potentially Visbilie View Cells" 
	// of the objects in front and back pvs

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

	interior->mTimeStamp = mTimeStamp ++;
	
	return interior;
}


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

	// remove the parents from the object pvss
	ObjectPvsMap::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
	ObjectPvsMap::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);
	}
}


bool VspTree::AddKdLeafToPvs(KdLeaf *leaf, 
							 ViewCell *vc, 
							 const float pdf, 
							 float &contribution)
{
	bool contri = false;

#if 1 // add kd intersecable to pvs
	KdIntersectable *kdObj = mOspTree->GetOrCreateKdIntersectable(leaf);
	
	if (vc->AddPvsSample(kdObj, pdf, contribution))
	{
		return true;
	}

#else // add all objects of kd node

	contribution = 0;

	ObjectContainer::const_iterator it, it_end = leaf->mObjects.end();

	for (it = leaf->mObjects.begin(); it != it_end; ++ it)
	{
		Intersectable *object = *it;						

		float newcontri;
						
		if (vc->AddPvsSample(object, pdf, newcontri))
		{
			contri = true;
		}

		//pdf += newPdf;
		newContri += contribution;
	}

#endif

	return contri;
}

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;

		Intersectable *obj = ray->mTerminationObject;

		if (obj) 
		{
			if (mStoreKdPvs)
			{
				// potentially visible kd cells
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
				AddKdLeafToPvs(leaf, vc, ray->mPdf, contribution);
			}
			else
			{
				if (vc->AddPvsSample(obj, ray->mPdf, contribution))
				{
					madeContrib = true;
				}
			}

			sc += contribution;
		}

		obj = ray->mOriginObject;

		if (obj) 
		{
			// potentially visible kd cells
			if (!mStoreKdPvs) // matt Todo: remove storekdpvs!!
			{
				if (vc->AddPvsSample(obj, ray->mPdf, contribution))
				{
					madeContrib = true;
				}
			}
			else
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);	
				AddKdLeafToPvs(leaf, vc, ray->mPdf, contribution);
			}
			
			

			sc += contribution;
		}

		if (madeContrib)
		{
			++ contributingSamples;
		}

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


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

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

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

	mLocalSplitCandidates->reserve(requestedSize);


	float pos;

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

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

		pos = (*rit).ExtrapTermination(axis);

		mLocalSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin, 
									pos, (*rit).mRay));
	}

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


int VspTree::GetPvsContribution(Intersectable *object) const
{
    int pvsContri = 0;

	KdPvsMap::const_iterator kit, kit_end = object->mKdPvs.mEntries.end();

	Intersectable::NewMail();

	// Search kd leaves this object is attached to
	for (kit = object->mKdPvs.mEntries.begin(); kit != kit_end; ++ kit)
	{
		KdNode *l = (*kit).first;

		// new object found during sweep 
		// => increase pvs contribution of this kd node
		if (!l->Mailed())
		{
			l->Mail();
			++ pvsContri;
		}
	}

	return pvsContri;
}


int VspTree::PrepareHeuristics(KdLeaf *leaf)
{	
	int 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());
		//Debug << "adding " << (int)leaf->mObjects.size() << " " << leaf->mMultipleObjects.size() << endl;
	}
	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;
}


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

	int 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;
		Intersectable *oObject = ray->mOriginObject;

		if (oObject)
		{
			if (!mUseKdPvsForHeuristics)
			{
				if (!oObject->Mailed())
				{
					oObject->Mail();
					oObject->mCounter = 1;

					++ pvsSize;
				}
				else
				{
					++ oObject->mCounter;
				}
			}
			else
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);
				pvsSize += PrepareHeuristics(leaf);	
			}	
		}

		Intersectable *tObject = (*ri).mRay->mTerminationObject;

		if (tObject)
		{
			if (!mUseKdPvsForHeuristics)
			{
				if (!tObject->Mailed())
				{
					tObject->Mail();
					tObject->mCounter = 1;
					++ pvsSize;
				}
				else
				{
					++ tObject->mCounter;
				}
			}
			else
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
				pvsSize += PrepareHeuristics(leaf);	
			}	
		}
	}

	return pvsSize;
}


void VspTree::RemoveContriFromPvs(KdLeaf *leaf, int &pvs) const
{
	// leaf falls out of right pvs
	if (-- leaf->mCounter == 0)
	{
		pvs -= ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
	}

	//-- handle objects which are in several kd leaves separately 
	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;
		}
	}
}


void VspTree::AddContriToPvs(KdLeaf *leaf, int &pvs) const
{
	if (leaf->Mailed())
		return;
	
	leaf->Mail();

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

	//-- handle objects of several kd leaves separately
	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;
		}
	}					
}


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

	Intersectable *oObject = ray->mOriginObject;
	Intersectable *tObject = ray->mTerminationObject;

	if (oObject)
	{	
		if (!mUseKdPvsForHeuristics)
		{
			if (ci.type == SortableEntry::ERayMin)
			{
				if (!oObject->Mailed())
				{
					oObject->Mail();
					++ pvsLeft;
				}
			}
			else if (ci.type == SortableEntry::ERayMax)
			{
				if (-- oObject->mCounter == 0)
					-- pvsRight;
			}
		}
		else // per kd node
		{
			KdLeaf *node = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);

			// add contributions of the kd nodes
			if (ci.type == SortableEntry::ERayMin)
			{
				AddContriToPvs(node, pvsLeft);
			}
			else if (ci.type == SortableEntry::ERayMax)
			{
				RemoveContriFromPvs(node, pvsRight);
			}
		}
	}
	
	if (tObject)
	{	
		if (!mUseKdPvsForHeuristics)
		{
			if (ci.type == SortableEntry::ERayMin)
			{
				if (!tObject->Mailed())
				{
					tObject->Mail();
					++ pvsLeft;
				}
			}
			else if (ci.type == SortableEntry::ERayMax)
			{
				if (-- tObject->mCounter == 0)
					-- pvsRight;
			}
		}
		else // per kd node
		{
			KdLeaf *node = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);

			if (ci.type == SortableEntry::ERayMin)
			{
				AddContriToPvs(node, pvsLeft);
			}
			else if (ci.type == SortableEntry::ERayMax)
			{
				RemoveContriFromPvs(node, pvsRight);
			}
		}
	}
}


float VspTree::EvalLocalCostHeuristics(const VspTraversalData &tData,
									   const AxisAlignedBox3 &box,
									   const int axis,
									   float &position)
{
	RayInfoContainer usedRays;

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

	int pvsSize = tData.mPvs;

	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;

	SortSplitCandidates(usedRays, axis, minBand, maxBand);

	// prepare the sweep
	// (note: returns pvs size, so there would be no need 
	// to give pvs size as argument)
	pvsSize = 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

	int pvsl = 0;
	int pvsr = pvsSize;

	int pvsBack = pvsl;
	int pvsFront = pvsr;

	float sum = (float)pvsSize * 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();


	//-- traverse through visibility events

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

	for (ci = mLocalSplitCandidates->begin(); ci != ci_end; ++ ci)
	{
		EvalPvsIncr(*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);
			//Debug  << "pos=" << (*ci).value << "\t newpos=" << currentPos << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << "\t cost= " << sum << endl;

			if (sum < minSum)
			{
				splitPlaneFound = true;

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

	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();

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

	const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
	
	const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

	if (splitPlaneFound)
	{
		ratio = newRenderCost / oldRenderCost;
	}
	
	const float volRatio = tData.mBoundingBox.GetVolume() / (sizeBox * mBoundingBox.GetVolume());

/*	
//if (axis != 1)
	//Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
	//      <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;

Debug << "\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;
*/
	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;

	// if we use some kind of specialised fixed axis
    const bool useSpecialAxis = 
		mOnlyDrivingAxis || mCirculatingAxis;
	//Debug << "data: " << tData.mBoundingBox << " pvs " << tData.mPvs << endl;
	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]);			
			}
			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;
	plane.mPosition = nPosition[bestAxis]; // split plane position

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

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


float VspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
									  const VspTraversalData &data,
									  float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;
	float totalPvs = 0;

	// probability that view point lies in back / front node
	float pOverall = data.mProbability;
	float pFront = 0;
	float pBack = 0;

	const float viewSpaceVol = mBoundingBox.GetVolume();

	// create unique ids for pvs heuristics
	Intersectable::NewMail(3);
	KdLeaf::NewMail(3);

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

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

		float t;
		VssRay *ray = rayInf.mRay;

		const int cf = 
			rayInf.ComputeRayIntersection(candidatePlane.mAxis, 
										  candidatePlane.mPosition, t);

		if (!mUseKdPvsForHeuristics)
		{
			// find front and back pvs for origing and termination object
			UpdateObjPvsContri(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
			UpdateObjPvsContri(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
		}
		else
		{
			if (ray->mTerminationObject)
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
				UpdateKdLeafPvsContri(leaf, cf, pvsFront, pvsBack, totalPvs);
			}

			if (ray->mOriginObject)
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);
				UpdateKdLeafPvsContri(leaf, cf, pvsFront, pvsBack, totalPvs);
			}
		}
	}


	AxisAlignedBox3 frontBox;
	AxisAlignedBox3 backBox;

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

	pFront = frontBox.GetVolume();
	pBack = pOverall - pFront;
		

	//-- pvs rendering heuristics
	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();

	//-- only render cost heuristics or combined with standard deviation
	const float penaltyOld = EvalPvsPenalty((int)totalPvs, lowerPvsLimit, upperPvsLimit);
    const float penaltyFront = EvalPvsPenalty((int)pvsFront, lowerPvsLimit, upperPvsLimit);
	const float penaltyBack = EvalPvsPenalty((int)pvsBack, lowerPvsLimit, upperPvsLimit);
			
	const float oldRenderCost = pOverall * penaltyOld;
	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

	const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;
	/*
	Debug << "\n==== eval 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;
*/
	return renderCostDecrease;
}


float VspTree::EvalLocalSplitCost(const VspTraversalData &data,
								  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);

	const int pvsSize = data.mPvs;

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

	// this is the main ray classification loop!
	for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
	{
		// determine the side of this ray with respect to the plane
		float t;
		const int side = (*rit).ComputeRayIntersection(axis, position, t);
	
		UpdateObjPvsContri((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
		UpdateObjPvsContri((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
	}

	//-- evaluate cost heuristics

	float pOverall;
	
	pOverall = data.mProbability;

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

	return  (mCtDivCi + newCost) / oldCost;
}


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

	//const float renderCost = mViewCellsManager->SimpleRay &raynderCost(obj);
	const int renderCost = 1;

	// 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::UpdateKdLeafPvsContri(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;
	}

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

	for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
	{	
		UpdateObjPvsContri(*oit, cf, frontPvs, backPvs, totalPvs);
    }	
	
	// 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 (!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().CountObjectsInPvs() <= 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 = data.mPvs;
	}

	mVspStats.pvs += data.mPvs;

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

	// accumulate rays to compute rays /  leaf
	mVspStats.accumRays += (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 _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().CountObjectsInPvs() << "), "
		  << "#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
				GetOrCreateOutOfBoundsCell()->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()
{
	mVspStats.invalidLeaves = 0;

	stack<VspNode *> nodeStack;

	if (!mRoot)
		return;

	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
{
	stack<VspNode *> nodeStack;

	if (!root)
		return;

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


void VspTree::CollectPvs(const RayInfoContainer &rays, 
						 ObjectContainer &objects) const
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	Intersectable::NewMail();

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

		Intersectable *object;
		object = ray->mOriginObject;

        if (object)
		{
			if (!object->Mailed())
			{
				object->Mail();
				objects.push_back(object);
			}
		}

		object = ray->mTerminationObject;

		if (object)
		{
			if (!object->Mailed())
			{
				object->Mail();
				objects.push_back(object);
			}
		}
	}
}


int VspTree::EvalPvsSize(const RayInfoContainer &rays) const
{
	int pvsSize = 0;
	Intersectable::NewMail();

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

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

			if (ray->mOriginObject)
			{
				if (!ray->mOriginObject->Mailed())
				{
					ray->mOriginObject->Mail();
					++ pvsSize;
				}
			}
			if (ray->mTerminationObject)
			{
				if (!ray->mTerminationObject->Mailed())
				{
					ray->mTerminationObject->Mail();
					++ pvsSize;
				}
			}
		}
	}
	else // compute pvs from kd nodes
	{
		KdNode::NewMail();

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

			if (ray->mTerminationObject)
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);

				if (!leaf->Mailed())
				{
					leaf->Mail();
					pvsSize += (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();
							++ pvsSize;
						}
					}
				}
			}

			if (ray->mOriginObject)
			{
				KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);

				if (!leaf->Mailed())
				{
					leaf->Mail();
					pvsSize += (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();
							++ pvsSize;
						}
					}
				}
			}
		}
	}

	return pvsSize;
}


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


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

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

	stack<LineTraversalData> tStack;

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

	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 *vc = leaf->GetViewCell();

			// don't have to mail because each view cell belongs to exactly one leaf
			//if (!vc->Mailed())
			//{
			//	vc->Mail();
				viewcells.push_back(vc);
				++ hits;
			//}
#if 0
			leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
#endif
			// 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;
			}
#if 0
			leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
#endif
			// 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()) 
		{
			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->GetBack());
			else
				nodeStack.push(interior->GetFront());
		}
	}
  
	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();

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

		const AxisAlignedBox3 bbox = GetBoundingBox(node);

		if (bbox.Includes(box))
		{
			// node geometry is contained in box
			CollectViewCells(node, true, viewCells, true);
		}
		else if (Overlap(bbox, box))
		{
			if (node->IsLeaf())
			{
				BspLeaf *leaf = dynamic_cast<BspLeaf *>(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::CollectDirtyCandidates(VspSplitCandidate *sc, 
									 vector<SplitCandidate *> &dirtyList)
{
	VspTraversalData &tData = sc->mParentData;
	VspLeaf *node = tData.mNode;
	
	KdLeaf::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;
	
		KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);
	
		if (!leaf->Mailed())
		{
			leaf->Mail();
			dirtyList.push_back(leaf->mSubdivisionCandidate);
		}
		
		leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
	
		if (!leaf->Mailed())
		{
			leaf->Mail();
			dirtyList.push_back(leaf->mSubdivisionCandidate);
		}
	}
}


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)
{
	/*mViewCellsManager->ComputeSampleContribution(ray, false, true);
	viewCells = ray.mViewCells;
	ray.mViewCells.clear();
	*/
	static Ray hray;
	hray.Init(ray);
	//hray.mFlags |= Ray::CULL_BACKFACES;
	//Ray hray(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);

	// if no precomputation of view cells
	CastLineSegment(origin, termination, viewCells);
}



/*******************************************************************/
/*                  class OspTree implementation                   */
/*******************************************************************/


OspTree::OspTree():
mRoot(NULL),
mTimeStamp(1),
mCopyFromKdTree(false)
{
	ReadEnvironment();
	mSplitCandidates = new vector<SortableEntry>;
}


OspTree::OspTree(const KdTree &kdTree)
{
	ReadEnvironment();

	// copy values from kd tree
	mCopyFromKdTree = true;
	mBoundingBox = kdTree.GetBox();
	mRoot = kdTree.GetRoot();

	mSplitCandidates = new vector<SortableEntry>;
}


OspTree::~OspTree()
{
	// delete kd intersectables
	KdIntersectableMap::iterator it, it_end = mKdIntersectables.end();

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

	// if not using kd tree root, delete tree (otherwise mKdTree has to do the job)
	if (!mCopyFromKdTree)
		DEL_PTR(mRoot);

	DEL_PTR(mSplitCandidates);
	mSubdivisionStats.close();
}


void OspTree::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", mTermMinProbability);
	
	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 << "******* OSP tree options ******** " << endl;

    Debug << "max depth: " << mTermMaxDepth << endl;
	Debug << "min probabiliy: " << mTermMinProbability << endl;
	Debug << "min objects: " << mTermMinObjects << endl;
	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	Debug << "miss tolerance: " << mTermMissTolerance << endl;
	Debug << "max leaves: " << mTermMaxLeaves << endl;

	Debug << "randomize: " << randomize << endl;

	Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
	Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
	Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
	Debug << "max memory: " << mMaxMemory << endl;
	Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
	Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
	
	Debug << "split borders: " << mSplitBorder << endl;
	Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
	Debug << endl;
}


void OspTree::SplitObjects(KdLeaf *parent,
						   const AxisAlignedPlane& splitPlane,
						   const ObjectContainer &objects,
						   ObjectContainer &front,
						   ObjectContainer &back)
{
	ObjectContainer::const_iterator oit, oit_end = objects.end();

    for (oit = objects.begin(); oit != oit_end; ++ oit) 
	{
		Intersectable *object = *oit;
		
		// initialise leaf references of objects
		if (parent->IsRoot())
		{
			object->mReferences = 1;
		}

		// remove parent ref
		-- object->mReferences;

		// determine the side of this ray with respect to the plane
		const AxisAlignedBox3 box = object->GetBox();

		if (box.Max(splitPlane.mAxis) >= splitPlane.mPosition)
		{
            front.push_back(object);
			++ object->mReferences;
		}

		if (box.Min(splitPlane.mAxis) < splitPlane.mPosition)
		{
			back.push_back(object);
			++ object->mReferences;
		}
	}

	mOspStats.objectRefs -= (int)objects.size();
	mOspStats.objectRefs += (int)(back.size() + front.size());
}


KdInterior *OspTree::SubdivideNode(
								   const AxisAlignedPlane &splitPlane,
								   const OspTraversalData &tData,
								   OspTraversalData &frontData,
								   OspTraversalData &backData
								   )
{
	KdLeaf *leaf = tData.mNode;
	
	// two new leaves
    mOspStats.nodes += 2;

	// add the new nodes to the tree
	KdInterior *node = new KdInterior(leaf->mParent);

	const int axis = splitPlane.mAxis;
	const float position = splitPlane.mPosition;

	node->mAxis = axis;
	node->mPosition = position;
	node->mBox = tData.mBoundingBox;


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

	// bounding boxes: split front and back node geometry
	tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition, 
		frontData.mBoundingBox, backData.mBoundingBox);

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

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

	//-- classify objects

	int objectsBack = 0;
	int objectsFront = 0;

    ObjectContainer::const_iterator mi, mi_end = leaf->mObjects.end();

	for ( mi = leaf->mObjects.begin(); mi != mi_end; ++ mi) 
	{
		// determine the side of this ray with respect to the plane
		const AxisAlignedBox3 box = (*mi)->GetBox();
		
		if (box.Max(axis) >= position) 
			++ objectsFront;
    
		if (box.Min(axis) < position)
			++ objectsBack;
	}


	// TODO matt: compute pvs
	//frontData.mPvs = objectsFront;
	//backData.mPvs = objectsBack;

	//CheckViewCellsPvs(leaf, *tData.mRays);
	RemoveParentViewCellsPvs(leaf, *tData.mRays);


	KdLeaf *back = new KdLeaf(node, objectsBack);
	KdLeaf *front = new KdLeaf(node, objectsFront);

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

	KdInterior *parent = leaf->mParent;

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

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

	SplitObjects(leaf, splitPlane, leaf->mObjects, front->mObjects, back->mObjects);
	
	FilterRays(front, *tData.mRays, *frontData.mRays);
	FilterRays(back, *tData.mRays, *backData.mRays);

	++ mOspStats.splits[splitPlane.mAxis];

	//-- eval view cells pvs
	
	// remove parent pvs update pvs of left and right leaves
	// Note: It is necessary to update first left and then right pvs.
	// We need to store the view cells seen by each object,
	// but also if the view cells are seen as part of two different
	// kd leaves, which is stored in the pdf component of the pvs entry. 
	// Because if an object is part of a view cell more than once,
	// it cannot be removed from the pvs by splitting a kd node where
	// the view cell sees only the other child of the node.
	// This is important during the subdivision

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

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


	ProcessMultipleRefs(front);
    ProcessMultipleRefs(back);

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


	// compute probability, i.e., volume of seen view cells
	frontData.mProbability = EvalViewCellsVolume(front, *frontData.mRays);
	backData.mProbability =  EvalViewCellsVolume(back, *backData.mRays);


	//delete leaf;
	return node;
}


KdNode *OspTree::Subdivide(SplitQueue &tQueue,
						   SplitCandidate *splitCandidate,
						   const bool globalCriteriaMet)
{
	OspSplitCandidate *sc = dynamic_cast<OspSplitCandidate *>(splitCandidate);
	OspTraversalData &tData = sc->mParentData;

	KdNode *newNode = tData.mNode;

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

		//-- continue subdivision
		
		// create new interior node and two leaf node
		const AxisAlignedPlane splitPlane = sc->mSplitPlane;
		
		newNode = SubdivideNode(splitPlane, 
								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
		
		OspSplitCandidate *frontCandidate = new OspSplitCandidate(tFrontData);
		OspSplitCandidate *backCandidate = new OspSplitCandidate(tBackData);

		EvalSplitCandidate(*frontCandidate);
		EvalSplitCandidate(*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)
		{
			KdLeaf *leaf = dynamic_cast<KdLeaf *>(newNode);

			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);
			}
		}
	}
	
	//-- cleanup
	tData.Clear();
	
	return newNode;
}


void OspTree::EvalSplitCandidate(OspSplitCandidate &splitCandidate)
{
	float frontProb;
	float backProb;

	// compute locally best split plane
	const bool success = 
		SelectSplitPlane(splitCandidate.mParentData, 
						 splitCandidate.mSplitPlane, 
						 frontProb, 
						 backProb);

	float oldRenderCost;

	// compute global decrease in render cost
	const float renderCostDecr = EvalRenderCostDecrease(splitCandidate.mSplitPlane, 
														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 OspTree::LocalTerminationCriteriaMet(const OspTraversalData &data) const
{
	// matt: TODO
	return ( 
		//(data.mNode->mObjects.size() < mTermMinObjects) ||
		//(data.mProbability <= mTermMinProbability) ||
		(data.mDepth >= mTermMaxDepth)
		 );
}


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


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

	/*if (data.mPvs > mOspStats.maxPvs)
	{
		mOspStats.maxPvs = data.mPvs;
	}

	mOspStats.pvs += data.mPvs;

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

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

#ifdef _DEBUG
	Debug << "BSP stats: "
		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
	  	  << "Prob: " << data.mProbability << " (min: " << mTermMinProbability << ")\n";
#endif
}


float OspTree::EvalLocalCostHeuristics(const OspTraversalData &tData,
									   const AxisAlignedBox3 &box,
									   const int axis,
									   float &position,
									   int &objectsFront,
									   int &objectsBack)
{
	RayInfoContainer usedRays;
	int mMaxTests = 500000; // HACK

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

	// 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 = box.Min(axis);
	const float maxBox = box.Max(axis);
	Debug << "min: " << minBox << " max: " << maxBox << endl;

	const float sizeBox = maxBox - minBox;

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

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

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

	float totalVol = 0;
	float voll = 0;
	float volr = totalVol;

	ViewCellContainer touchedViewCells;
	
	const float totalRenderCost = PrepareHeuristics(tData, touchedViewCells);
	float renderCost = totalRenderCost;

	/// this is kind of a reverse pvs
	const int totalPvs = (int)tData.mNode->mObjects.size();
	
	int pvsl = 0;
	int pvsr = totalPvs;

	float sum = (float)totalVol * sizeBox;


	Debug << "here82 render cost: " << renderCost / viewSpaceVol << endl;

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

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

	float volBack = voll;
	float volFront = volr;

	int pvsBack = pvsl;
	int pvsFront = pvsr;
	
	float minRenderCost = 1e20f;

	debugVol = 0;



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

	
	//-- traverse through events and find best split plane
	
	vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
	//minRenderCost = RandomValue(0,viewSpaceVol);

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

		EvalHeuristicsContribution(tData.mNode, 
								   *ci, 
								   renderCost,
								   touchedViewCells
								  );

		// Note: sufficient to compare size of bounding boxes of front and back side?
		if (((*ci).mPos >= minBand) && ((*ci).mPos <= maxBand))
		{
			float currentPos;
			//Debug << "here29 : " << minRenderCost / viewSpaceVol << endl;

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

			if (renderCost < minRenderCost)
			{
				splitPlaneFound = true;
				minRenderCost = renderCost;
				position = currentPos;
				Debug << "pos: " << position << endl;
			}
		}
	}
//splitPlaneFound = true;
	if (splitPlaneFound)
	{
		ratio = minRenderCost / totalRenderCost;
	}

	Debug << "\n§§§§ eval local cost §§§§" << endl
		  << "old rc: " << totalRenderCost / viewSpaceVol << " new rc: " << minRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << (totalRenderCost - minRenderCost) / viewSpaceVol  << endl;

	return ratio;
}


void OspTree::SortSplitCandidates(const OspTraversalData &tData,
								  const int axis, 
								  float minBand, 
								  float maxBand)

{
	mSplitCandidates->clear();

	RayInfoContainer *rays = tData.mRays;
	KdLeaf *leaf = tData.mNode;

	int requestedSize = 2 * (int)rays->size() + 2 * (int)leaf->mObjects.size();

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

	mSplitCandidates->reserve(requestedSize);

	float pos;

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

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

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

		const bool positive = ray->HasPosDir(axis);
	
		if (EndPointInsideNode(leaf, *ray, true))
		{
			pos = ray->mTermination[axis];

			mSplitCandidates->push_back(
				SortableEntry(SortableEntry::BOX_INTERSECT, 
							  pos, 
							  ray->mTerminationObject, 
							  ray)
							  );
		}

		// if hitpoint with object is inside this node
		if (EndPointInsideNode(leaf, *ray, false))
		{
			pos = ray->mOrigin[axis];
	
			mSplitCandidates->push_back(
				SortableEntry(SortableEntry::BOX_INTERSECT, 
							  pos, 
							  ray->mOriginObject, 
							  ray)
							  );
		}
	}


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

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

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

		mSplitCandidates->push_back(
			SortableEntry(SortableEntry::BOX_MIN,
						  box.Min(axis),
						  object,
						  NULL)
						  );
   
    
		mSplitCandidates->push_back(
			SortableEntry(SortableEntry::BOX_MAX,
						  box.Max(axis),
						  object,
						  NULL)
						  );
	}

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


const OspTreeStatistics &OspTree::GetStatistics() const
{
	return mOspStats;
}


void OspTree::PrepareHeuristics(const VssRay &ray, 
								ViewCellContainer &touchedViewCells)
{
	// collect view cells and set mail + counter
	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())
		{
			vc->Mail(); // set as belonging to front leaf
			vc->mCounter = 0;
			touchedViewCells.push_back(vc);
		}
		
		// view cell volume already added => just increase counter
		++ vc->mCounter;
	}
}


float OspTree::PrepareHeuristics(const OspTraversalData &tData,
								 ViewCellContainer &touchedViewCells)
{	
	float renderCost = 0;

	Intersectable::NewMail(3);
	ViewCell::NewMail(3);
	MailablePvsData::NewMail();

	KdLeaf *leaf = tData.mNode;

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

	for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
     	PrepareHeuristics(*ray, touchedViewCells);
	}

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

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		obj->Mail(); // set as belonging to front leaf
		
		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) // should never happen
				continue;

			if (!vdata->Mailed())
			{
				vdata->Mail();
				renderCost += vc->GetVolume();
			}
		}
	}

	Debug << "here32 " << touchedViewCells.size() << endl;
	return renderCost;
}


void OspTree::AddObjectContribution(KdLeaf *leaf,
						  Intersectable *obj,
						  ViewCellContainer &touchedViewCells,
						  float &renderCost)
{
	//Debug << "add contri to obj: " << obj->mMailbox - Intersectable::sMailId << endl;
	obj->Mail(2); // set as belonging to both leafs

	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();
	
	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;
		//Debug << "here30 vc mail " << vc->mMailbox - ViewCell::sMailId << endl;

		// if obj not previously associated with this view cell => increase render cost
		if (vc->Mailed(1) && !ViewCellHasMultipleReferences(obj, vc, false))
		{	
			//Debug << "add to rendercost: " << vc->GetVolume() / mVspTree->GetBoundingBox().GetVolume() << endl;
			renderCost += vc->GetVolume();
		}
	}
}


void OspTree::SubtractObjectContribution(KdLeaf *leaf,
								Intersectable * obj,
								ViewCellContainer &touchedViewCells,
								float &renderCost)
{
	//Debug << "subtract contri from obj: " << obj->mMailbox - Intersectable::sMailId << endl;
    obj->Mail(1); // set as belonging to back leaf
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();
    	
	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;

		//Debug << "here6 " << vc->mMailbox - ViewCell::sMailId << endl;

		// if obj was previously associated with this view cell but is not now 
		// => decrease render cost
		if (vc->Mailed() &&	!ViewCellHasMultipleReferences(obj, vc, false))
		{
			//Debug << "subtract from rendercost: " << vc->GetVolume() / mVspTree->GetBoundingBox().GetVolume() << endl;
			renderCost -= vc->GetVolume();
		}
	}
}


void OspTree::EvalRayContribution(KdLeaf *leaf, const VssRay &ray, float &renderCost)
{
	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)
	{
		EvalViewCellContribution(leaf, *vit, renderCost);
	}
}


void OspTree::EvalViewCellContribution(KdLeaf *leaf, 
									   ViewCell *viewCell,
									   float &renderCost)
{
	//Debug << "**************" << endl;
	//Debug << "vc contri: " << renderCost / mVspTree->GetBoundingBox().GetVolume() << " " << viewCell->mMailbox - ViewCell::sMailId << " counter: " << viewCell->mCounter << endl;

	const float vol = viewCell->GetVolume();

	if (viewCell->Mailed())
	{
		viewCell->Mail(2); // view cell can be seen from both nodes

        // we now see view cell from both nodes => add contri
		ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

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

			// was render cost already added?
			if (!ViewCellHasMultipleReferences(obj, viewCell, false) &&
				obj->Mailed(1))
			{
				//Debug << "obj mail 1!! " << vol / mVspTree->GetBoundingBox().GetVolume() << endl;
				renderCost += vol;
			}
		}
	}

	if (-- viewCell->mCounter == 0)
	{
		viewCell->Mail(1); // view cell can be seen from back node only

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

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

			// can render cost be be reduced?
			if (!ViewCellHasMultipleReferences(obj, viewCell, false) &&	obj->Mailed())
			{
				renderCost -= vol;
			}
		}
	}
	//Debug << "new rc: " << renderCost / mVspTree->GetBoundingBox().GetVolume() << " counter2: " << viewCell->mCounter << endl;
	//Debug << "**************"<<endl;
}


void OspTree::EvalHeuristicsContribution(KdLeaf *leaf,
										 const SortableEntry &ci,
										 float &renderCost,
										 ViewCellContainer &touchedViewCells)
{
	Intersectable *obj = ci.mObject;
	VssRay *ray = ci.mRay;

	// switch between different types of events
	switch (ci.mType) 
	{
		case SortableEntry::BOX_MIN:
			cout << "<";
			AddObjectContribution(leaf, ci.mObject, touchedViewCells, renderCost);	
			break;
			
		case SortableEntry::BOX_MAX:
			cout << ">";
			SubtractObjectContribution(leaf, ci.mObject, touchedViewCells, renderCost);
			break;

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

		default:
			cout << "should not come here" << endl;
			break;
	}

	//cout << "vol left: " << volLeft << " vol right " << volRight << endl;
}


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


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


float OspTree::SelectSplitPlane(const OspTraversalData &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;

	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)
			{
				//-- place split plane using heuristics
				int objectsFront, objectsBack;

				nCostRatio[axis] =
					EvalLocalCostHeuristics(tData,
									   tData.mBoundingBox,
									   axis,
									   nPosition[axis],
									   objectsFront,
									   objectsBack);			
			}
			/*	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

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

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

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


bool OspTree::EndPointInsideNode(KdLeaf *leaf, 
								 VssRay &ray, 
								 bool isTermination) const
{
	// test if the leaf where the hitpoint is located is the current leaf
	if (isTermination)
	{
		 return ray.mTerminationObject && (GetLeaf(ray.mTermination, ray.mTerminationNode) == leaf);
	}
	else
	{
		return false; // no origin object!
		return ray.mOriginObject && (GetLeaf(ray.mOrigin, ray.mOriginNode) == leaf);
	}
}


void OspTree::EvalSubdivisionStats(const SplitCandidate &sc)
{
	const float costDecr = sc.GetRenderCostDecrease();

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


void OspTree::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;
}


int OspTree::ClassifyRay(VssRay *ray, 
						 KdLeaf *leaf, 
						 const AxisAlignedPlane &plane) const
{
	const bool originInside = EndPointInsideNode(leaf, *ray, false);
    const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

	const bool originGt = (ray->mOrigin[plane.mAxis] >= plane.mPosition);
	const bool originLt = (ray->mOrigin[plane.mAxis] <= plane.mPosition);

	const bool terminationGt = (ray->mTermination[plane.mAxis] >= plane.mPosition);
	const bool terminationLt = (ray->mTermination[plane.mAxis] <= plane.mPosition);

	// add view cell volume to front volume
	const bool addToFront = 
		((originInside && originGt) || (terminationInside && terminationGt));

	// add view cell volume to back volume
	const bool addToBack = 
		((originInside && originLt/*!originGt*/) || (terminationInside && terminationLt/*terminationGt*/));

	// classify ray with respect to the child nodes the view cells contribute
	if (addToFront && addToBack)
		return 0;
	else if (addToBack)
		return -1;
	else if (addToFront)
		return 1;

	return -2;
}


int OspTree::ClassifyRays(const RayInfoContainer &rays, 
						  KdLeaf *leaf, 
						  const AxisAlignedPlane &plane,
						  RayInfoContainer &frontRays,
						  RayInfoContainer &backRays) const
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

		bool originGt = false;
		bool terminationGt = false;
		
		Intersectable *obj = ray->mOriginObject;

		if (0 && obj)
		{
			originGt = ray->mOrigin[plane.mAxis] >= plane.mPosition;
		}

		obj = ray->mTerminationObject;
		
		if (obj)
		{
			terminationGt = ray->mTermination[plane.mAxis] >= plane.mPosition;
		}

		// either begin or end point on front side
		const bool addToFront = originGt || terminationGt;

		// add view cell volume to back volume
		const bool addToBack = !originGt || !terminationGt;
	
		// classify ray with respect to the child nodes the view cells contribute
		if (addToFront)
			frontRays.push_back(*rit);
		if (addToBack)
			backRays.push_back(*rit);
	}

	return 0;
}


#if 0

float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
									  const OspTraversalData &tData,
									  float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;

	// probability that view point lies in back / front node
	float pOverall = 0;//data.mProbability; // todo matt§: value ok?
	float pFront = 0;
	float pBack = 0;
	float pFrontAndBack = 0;

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

	Intersectable::NewMail(3);
	
	KdLeaf *leaf = tData.mNode;
	const int totalPvs = (int)leaf->mObjects.size();
	
	RayInfoContainer frontRays, backRays;
	ClassifyRays(*tData.mRays, leaf, candidatePlane, frontRays, backRays);

	ViewCellContainer touchedViewCells;
	
	// 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);

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

	for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// add volume to volumes of left and / or right children
		// if one of the ray end points is inside
		const int classification = ClassifyRay(ray, leaf, candidatePlane);

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

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// if not previously mailed
			if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
			{
				touchedViewCells.push_back(vc);
			}

			// classify / mail the view cell
			MailViewCell(vc, classification);	
		}
	}

	// evaluate view cells volume contribution 
	// with respect to the mail box classification
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
	}

	//////////////////////////////////////////////
	//
	// evaluate contribution of objects which are part of other kd nodes
	// which are seen by the same view cell
	// These contributions cannot be reduced by splitting, because
	// the object would still be part of the pvs

	float additionalFrontRenderCost = 0;
	float additionalBackRenderCost = 0;
	
	MailablePvsData::NewMail();

	for (rit = tData.mRays->begin(); rit != tData.mRays->end(); ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

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

		// traverse through view cells
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *viewCell = *vit;

			// for a view cell:
			// if object is also in the kd pvs entries from other kd leaves,
			// render cost cannot be reduced for this view cell
			// => the render cost was falsly reduced, add them again

			if (ray->mTerminationObject)
			{	
				Intersectable *obj = ray->mTerminationObject;
				
				const bool renderCostWrong = 
					ViewCellHasMultipleReferences(obj, viewCell, true);
			
				// if there is already an entry of this object in the view cell pvs 
				if (renderCostWrong)
				{
					// view cell was counted only for front or back => correct tjos
					if (viewCell->Mailed(1) && obj->Mailed())
					{
						additionalFrontRenderCost += viewCell->GetVolume();
					}
					else if (viewCell->Mailed() && obj->Mailed(1))
					{
						additionalBackRenderCost += viewCell->GetVolume();
					}
				}
			}
	
			if (ray->mOriginObject)
			{
				Intersectable *obj = ray->mOriginObject;

				const bool renderCostWrong = 
					ViewCellHasMultipleReferences(obj, viewCell, true);
				// if there is already an entry of this object in the view cell pvs 
				if (renderCostWrong)
				{
					if (viewCell->Mailed(1) && obj->Mailed())
					{
						additionalFrontRenderCost += viewCell->GetVolume();
					}
					else if (viewCell->Mailed() && obj->Mailed(1))
					{
						additionalBackRenderCost += viewCell->GetVolume();
					}
				}
			}
		}
	}

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


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


	//-- pvs rendering heuristics

	const float oldRenderCost = pOverall * totalPvs;
	
	// sum up the terms: 
	// the view cells seeing
	// a) the left node are weighted by the #left node objects
	// b) the right node are weighted by the #right node objects
	// c) both nodes are weighted by the #parent node objects
	const float newRenderCost = pvsFront * pFront + pvsBack * pBack + totalPvs * pFrontAndBack
		+ additionalFrontRenderCost + additionalBackRenderCost;

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

	Debug << "\n(((( eval render cost decrease ))))" << endl
		  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl 
		  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol 
		  << " front and back p " << pFrontAndBack / viewSpaceVol << " p: " << tData.mProbability / viewSpaceVol << endl
		  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << renderCostDecrease << endl
		  << "additional front " << additionalFrontRenderCost / viewSpaceVol 
		  << " additional back " << additionalBackRenderCost  / viewSpaceVol << endl;

	if (oldRenderCost < newRenderCost * 0.99)
		Debug <<"error!!"<<endl;

	//if ((((pOverall - tData.mProbability) / viewSpaceVol) > 0.00001))
	//	Debug << "ERROR!!"<<endl;

	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

		
	return renderCostDecrease;
}

#else

float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
									  const OspTraversalData &tData,
									  float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;

	// probability that view point lies in back / front node
	float pOverall = 0;//data.mProbability; // todo matt§: value ok?
	float pFront = 0;
	float pBack = 0;
	float pFrontAndBack = 0;

	Intersectable::NewMail(3);
	KdLeaf *leaf = tData.mNode;

	const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();
	const int totalPvs = (int)leaf->mObjects.size();
	
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.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 = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		int pvsContri = 0;

		Intersectable *obj = *oit;
		const AxisAlignedBox3 box = obj->GetBox();

		// test if box falls in left / right child node
		if (box.Max(candidatePlane.mAxis) >= candidatePlane.mPosition)
		{
			++ pvsFront;
			obj->Mail();		
		}

		if (box.Min(candidatePlane.mAxis) < candidatePlane.mPosition)
		{
			++ pvsBack;

			if (obj->Mailed())
				obj->Mail(2);
			else 
				obj->Mail(1);
		}
	}

	
	ViewCellContainer touchedViewCells;

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

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
	//map<Intersectable *, ViewCellContainer> objectsMap;

	for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// add volume to volumes of left and / or right children
		// if one of the ray end points is inside
		const int classification = ClassifyRay(ray, leaf, candidatePlane);

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

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// if not previously mailed
			if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
			{
				touchedViewCells.push_back(vc);
			}

			// classify / mail the view cell
			MailViewCell(vc, classification);
		}
	}

	// evaluate view cells volume contribution 
	// with respect to the mail box classification
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
	}

	//////////////////////////////////////////////
	//
	// evaluate contribution of objects which are part of other kd nodes
	// which are seen by the same view cell
	// These contributions cannot be reduced by splitting, because
	// the object would still be part of the pvs

	float newRc = 0;
	float rc = 0;

	MailablePvsData::NewMail();
	//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();

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

			if (!vdata)
			{
				Debug << "here86 error!!"<<endl;
				continue;
			}

			if (!vdata->Mailed())
			{
				vdata->Mail();
				rc += vc->GetVolume();

				if ((vdata->mSumPdf > 1.5) || 
					vc->Mailed(2) || 
					obj->Mailed(2) ||
					(obj->Mailed() && vc->Mailed()) ||
					(obj->Mailed(1) && vc->Mailed(1))
					)
				{
					newRc += vc->GetVolume();
				}
			}
		}
	}


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


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


	//-- pvs rendering heuristics

	const float oldRenderCost = rc;//pOverall * totalPvs;
	
	// sum up the terms: 
	// the view cells seeing
	// a) the left node are weighted by the #left node objects
	// b) the right node are weighted by the #right node objects
	// c) both nodes are weighted by the #parent node objects
	const float newRenderCost = newRc//pvsFront * pFront + pvsBack * pBack + totalPvs * pFrontAndBack
;//		+ additionalFrontRenderCost + additionalBackRenderCost;

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

	Debug << "\n(((( eval render cost decrease ))))" << endl
		  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl 
		  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol 
		  << " front and back p " << pFrontAndBack / viewSpaceVol << " p: " << tData.mProbability / viewSpaceVol << endl
		  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << renderCostDecrease << endl;
		  //<< "additional front " << additionalFrontRenderCost / viewSpaceVol 
		  //<< " additional back " << additionalBackRenderCost  / viewSpaceVol << endl;

	if (oldRenderCost < newRenderCost * 0.99)
		Debug <<"error!!"<<endl;

	//if ((((pOverall - tData.mProbability) / viewSpaceVol) > 0.00001))
	//	Debug << "ERROR!!"<<endl;

	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

		
	return renderCostDecrease;
}
#endif


void OspTree::ComputeBoundingBox(const ObjectContainer &objects,
								 AxisAlignedBox3 *forcedBoundingBox) 
{
	if (forcedBoundingBox)
	{
		mBoundingBox = *forcedBoundingBox;
	}
	else // compute vsp tree bounding box
	{
		mBoundingBox.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());
		}
	}
}


void OspTree::ProcessMultipleRefs(KdLeaf *leaf) const
{
	// find objects from multiple kd-leaves
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *object = *oit;
		
		if (object->mReferences > 1)
		{
			leaf->mMultipleObjects.push_back(object);
		}
	}
}


void OspTree::CollectLeaves(vector<KdLeaf *> &leaves) const
{
	stack<KdNode *> nodeStack;
	nodeStack.push(mRoot);

	while (!nodeStack.empty()) 
	{
		KdNode *node = nodeStack.top();
		nodeStack.pop();
		if (node->IsLeaf()) 
		{
			KdLeaf *leaf = (KdLeaf *)node;
			leaves.push_back(leaf);
		}
		else 
		{
			KdInterior *interior = (KdInterior *)node;
			nodeStack.push(interior->mBack);
			nodeStack.push(interior->mFront);
		}
	}
}


AxisAlignedBox3 OspTree::GetBoundingBox(KdNode *node) const
{
	if (!node->mParent)
		return mBoundingBox;

	if (!node->IsLeaf())
	{
		return (dynamic_cast<KdInterior *>(node))->mBox;		
	}

	KdInterior *parent = dynamic_cast<KdInterior *>(node->mParent);

	AxisAlignedBox3 box(parent->mBox);

	if (parent->mFront == node)
		box.SetMin(parent->mAxis, parent->mPosition);
    else
		box.SetMax(parent->mAxis, parent->mPosition);

	return box;
}


void OspTree::CollectViewCells(KdLeaf *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 OspTree::CollectDirtyCandidates(OspSplitCandidate *sc, 
									 vector<SplitCandidate *> &dirtyList)
{
	OspTraversalData &tData = sc->mParentData;
	KdLeaf *node = tData.mNode;
	
	ViewCell::NewMail();
	ViewCellContainer viewCells;
	ViewCellContainer tmpViewCells;

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

	// find all view cells associated with the rays, add them to view cells
	for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
		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->GetSplitCandidate());
	}
}


KdNode *OspTree::GetRoot() const
{
	return mRoot;
}


KdLeaf *OspTree::GetLeaf(const Vector3 &pt, KdNode *node) const
{
	// start from root of tree
	if (node == NULL)
	{
		node = mRoot;
	}

	stack<KdNode *> nodeStack;
	nodeStack.push(node);
  
	KdLeaf *leaf = NULL;
  
	while (!nodeStack.empty())  
	{
		KdNode *node = nodeStack.top();
		nodeStack.pop();
	
		if (node->IsLeaf()) 
		{
			leaf = dynamic_cast<KdLeaf *>(node);
		} 
		else 	
		{	
			// find point
			KdInterior *interior = dynamic_cast<KdInterior *>(node);
        
			if (interior->mPosition > pt[interior->mAxis])
			{
				nodeStack.push(interior->mBack);
			}
			else
			{
				nodeStack.push(interior->mFront);
			}
		}
	}
  
	return leaf;
}


void OspTree::PreprocessRays(KdLeaf *root,
							 const VssRayContainer &sampleRays,
							 RayInfoContainer &rays)
{
	VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
	RayInfoContainer clippedRays;

	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;

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

			// HACK: reset nodes for the case we have used object kd tree for
			// tree building.
			ray->mOriginNode = NULL;
			ray->mTerminationNode = NULL;
		}
	}

	FilterRays(root, clippedRays, rays);

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



int OspTree::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;
}


int OspTree::FilterRays(KdLeaf *leaf, 
						const RayInfoContainer &rays,
						RayInfoContainer &filteredRays)
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

		// test if intersection point with one of the objects is inside this node
		const bool originInside = EndPointInsideNode(leaf, *ray, false);
        const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

		if (originInside || terminationInside)
		{
			filteredRays.push_back(ray);
		}
	}

	return 0;
}
					 

int OspTree::CheckViewCellsPvs(KdLeaf *leaf, 
							   const RayInfoContainer &rays) const
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	Intersectable::NewMail();
	ObjectContainer touchedObjects;
	

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

		if (ray->mTerminationObject)
		{
			Intersectable *obj = ray->mTerminationObject;

			if (!obj->Mailed())
			{
				obj->Mail();
				touchedObjects.push_back(obj);
			}
		}

		if (ray->mOriginObject)
		{
			Intersectable *obj = ray->mOriginObject;

			if (!obj->Mailed())
			{
				obj->Mail();
				touchedObjects.push_back(obj);
			}
		}
	}
	Debug << "here65 " << touchedObjects.size() << endl;
	ObjectContainer::const_iterator it, it_end = touchedObjects.end();
	for (it = touchedObjects.begin(); it != it_end; ++ it)
	{
		Debug << "\nhere94 obj: " << (*it) << " size: " << (*it)->mViewCellPvs.GetSize() << endl << endl;
		ViewCellPvsMap::const_iterator mit, mit_end = (*it)->mViewCellPvs.mEntries.end();

		for (mit = (*it)->mViewCellPvs.mEntries.begin(); mit != mit_end; ++ mit)
		{
			Debug << "newsumpdf: " << (*mit).second.mSumPdf << endl;
		}
	}

	return 0;
}


KdIntersectable *OspTree::GetOrCreateKdIntersectable(KdNode *node)
{
	// search nodes
	std::map<KdNode *, KdIntersectable *>::
		const_iterator it = mKdIntersectables.find(node);

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

	// not in map => create new entry
	KdIntersectable *kdObj = new KdIntersectable(node);
	mKdIntersectables[node] = kdObj;

	return kdObj;
}


/*
void OspTree::AddViewCellVolume(ViewCell *vc,
								const int cf,
								float &frontVol,
								float &backVol,
								float &totalVol) const
{
	//const float renderCost = mViewCellsManager->SimpleRay &raynderCost(obj);
	const float vol = vc->GetVolume();

	// view cell not found yet => new
	if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
	{
		totalVol += vol;
	}

	if (cf >= 0) // front volume
	{
		if (!vc->Mailed() && !vc->Mailed(2))
		{
			frontVol += vol;
		
			// already in back volume => in both volumes
			if (vc->Mailed(1))
				vc->Mail(2);
			else
				vc->Mail();
		}
	}

	if (cf <= 0) // back volume
	{
		if (!vc->Mailed(1) && !vc->Mailed(2))
		{
			backVol += vol;
		
			// already in front volume => in both volume
			if (vc->Mailed())
				vc->Mail(2);
			else
				vc->Mail(1);
		}
	}
}
*/


void OspTree::MailViewCell(ViewCell *vc, const int cf) const
{
	if (vc->Mailed(2)) // already classified
		return;

	if (cf >= 0) // front volume
	{
		// already in back volume => in both volumes
		if (vc->Mailed(1))
		{
			vc->Mail(2);
		}
		else
		{
			vc->Mail();
		}
	}

	if (cf <= 0) // back volume
	{
		// already in front volume => in both volume
		if (vc->Mailed())
		{
			vc->Mail(2);
		}
		else
		{
			vc->Mail(1);
		}
	}
}


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


float OspTree::EvalViewCellsVolume(KdLeaf *leaf, 
								   const RayInfoContainer &rays) const
{
	float vol = 0;
	ViewCell::NewMail();

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

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

		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())
			{
				vc->Mail();
				vol += vc->GetVolume();
			}
		}		
	}

	return vol;
}


bool OspTree::AddViewCellToObjectPvs(Intersectable *obj, 
									 ViewCell *vc, 
									 float &contribution,
									 bool onlyMailed) const
{	
	contribution = 0; // todo

	MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

	if (!vdata)
	{
		vdata = obj->mViewCellPvs.AddSample2(vc, 1);
	}
	else if (!onlyMailed || !vdata->Mailed())
	{
		obj->mViewCellPvs.AddSample(vc, 1);
	}
	
	vdata->Mail();

	return true;
}


void OspTree::CollectTouchedViewCells(const RayInfoContainer &rays, 
									  ViewCellContainer &touchedViewCells) const
{
	ViewCell::NewMail();
	
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

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

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			if (!vc->Mailed())
			{
				vc->Mail();
				touchedViewCells.push_back(vc);
			}
		}
	}
}


int OspTree::UpdateViewCellsPvs(KdLeaf *leaf, 
								const RayInfoContainer &rays) const

{
	MailablePvsData::NewMail();
	
	/*RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

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

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			Intersectable *obj = ray->mTerminationObject;
			if (obj)
			{
				float contri;
				AddViewCellToObjectPvs(obj, vc, contri, true);
			}

			obj = ray->mOriginObject;
			if (obj)
			{
				float contri;
				AddViewCellToObjectPvs(obj, vc, contri, true);
			}
		}
	}*/
	
	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 OspTree::RemoveParentViewCellsPvs(KdLeaf *leaf, 
									  const RayInfoContainer &rays
									  ) const

{
	MailablePvsData::NewMail();

	/*RayInfoContainer::const_iterator rit, rit_end = rays.end();
	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// test if intersection point with one of the objects is inside this node
		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

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

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			Intersectable *obj = ray->mTerminationObject;

			if (obj)
			{
				MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

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

			obj = ray->mOriginObject;

			if (obj)
			{
				MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

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

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


float OspTree::EvalRenderCost(const VssRayContainer &myrays)
{
	float rc = 0;
	float prop = mVspTree->GetBoundingBox().GetVolume();	

	KdLeafContainer leaves;
	CollectLeaves(leaves);

	ViewCellContainer vcleaves;
	mVspTree->CollectViewCells(vcleaves, false);
	

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

	for (vit = vcleaves.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;
		vc->GetPvs().Clear();
	}

	KdLeafContainer::const_iterator kit, kit_end = leaves.end();

	MailablePvsData::NewMail();

	for (kit = leaves.begin(); kit != kit_end; ++ kit)
	{
		KdLeaf *leaf = *kit;
		float newCost = 0;
		float vol = 0;

		ViewCell::NewMail();
		VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
		ViewCellContainer touchedViewCells;

		for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
		{
			VssRay *ray = *rit;
			
			// test if intersection point with one of the objects is inside this node
			const bool originInside = EndPointInsideNode(leaf, *ray, false);
			const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

			if (originInside || terminationInside)
			{
				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 not previously mailed
					if (!vc->Mailed())
					{
						vc->Mail();
						touchedViewCells.push_back(vc);
					}

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

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

				if (!vdata)
				{
					vc->GetPvs().AddSample(obj, 1);
					newCost += vc->GetVolume();
				}

/*				MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

				if (!vdata || !vdata->Mailed())
				{
					newCost += vc->GetVolume();
					if (!vdata)
						vdata = obj->mViewCellPvs.AddSample2(vc, 1);
					vdata->Mail();
				}*/
			}
		}

		rc += newCost;
	}

	return rc / prop;
}


float OspTree::EvalLeafCost(const OspTraversalData &tData)
{
	KdLeaf *leaf = tData.mNode;

	float rc = 0;
	//float prop = mVspTree->GetBoundingBox().GetVolume();	
	//float vol = 0;
	//ViewCell::NewMail();
	
	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
	ViewCellContainer touchedViewCells;

	for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
			
		// test if intersection point with one of the objects is inside this node

		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 not previously mailed
			if (!vc->Mailed())
			{
				vc->Mail();
				touchedViewCells.push_back(vc);
			}
		}
	}

	// clear pvs of involved view cells
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;
		vc->GetPvs().Clear();
	}
	
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	//Debug << "here53 " << touchedViewCells.size() << endl;
	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

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

			if (!vdata)
			{
				vc->GetPvs().AddSample(obj, 1);
				rc += vc->GetVolume();
				//prop += vc->GetVolume();
			}
		}
	}
	
	return rc;
}


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

	return true;
}


void OspTree::ExportNode(KdNode *node, OUT_STREAM &stream)
{
	if (node->IsLeaf())
	{
		KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);

		stream << "<Leaf ";
		stream << "objects=\"";
		
		//-- export objects in kd leaves
		//if (mExportObjects) ExportObjects(leaf, stream);
		
		stream << "\" />" << endl;
		//stream << " </Leaf>" << endl;
	}
	else
	{	
		KdInterior *interior = dynamic_cast<KdInterior *>(node);
	
		stream << "<Interior plane=\"" << interior->mPosition << " " 
			   << interior->mAxis << "\">" << endl;

		ExportNode(interior->mBack, stream);
		ExportNode(interior->mFront, stream);

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


struct KdObjectsTraversalData
{
	KdNode *node;
	ObjectContainer *objects;
};


void OspTree::InsertObjects(KdNode *node, const ObjectContainer &objects)
{
	stack<KdObjectsTraversalData> tStack;

	while (!tStack.empty())
	{
		KdObjectsTraversalData tData = tStack.top();
        tStack.pop();

		KdNode *node = tData.node;
		
		if (node->IsLeaf())
		{
			KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);

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

			for (oit = tData.objects->begin(); oit != oit_end; ++ oit)
			{
				leaf->mObjects.push_back(*oit);
			}
		}
		else // interior
		{
			KdObjectsTraversalData frontData, backData;
			KdInterior *interior = dynamic_cast<KdInterior *>(node);

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

			ObjectContainer::const_iterator oit, oit_end = tData.objects->end();
			
		    for (oit = tData.objects->begin(); oit != oit_end; ++ oit) 
			{
				Intersectable *object = *oit;
		
				// determine the side of this ray with respect to the plane
				const AxisAlignedBox3 box = object->GetBox();

				if (box.Max(interior->mAxis) >= interior->mPosition)
				{
					frontData.objects->push_back(object);
				}

				if (box.Min(interior->mAxis) < interior->mPosition)
				{
					backData.objects->push_back(object);
				}
			}

			tStack.push(backData);
			tStack.push(frontData);
		}

		DEL_PTR(tData.objects);
	}
}


/*******************************************************************/
/*              class HierarchyManager implementation              */
/*******************************************************************/


HierarchyManager::HierarchyManager(VspTree &vspTree, OspTree &ospTree):
mVspTree(vspTree), mOspTree(ospTree)
{
	// cross references
	mVspTree.mOspTree = &ospTree;
	mOspTree.mVspTree = &vspTree;

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


SplitCandidate *HierarchyManager::NextSplitCandidate()
{
	SplitCandidate *splitCandidate = mTQueue.Top();
	mTQueue.Pop();

	return splitCandidate;
}


VspTree::VspSplitCandidate *HierarchyManager::PrepareVsp(const VssRayContainer &sampleRays,
											 AxisAlignedBox3 *forcedViewSpace,
											 RayInfoContainer &rays)
{
	// store pointer to this tree
	VspTree::VspSplitCandidate::sVspTree = &mVspTree;
	mVspTree.mVspStats.nodes = 1;

	// compute view space bounding box
	mVspTree.ComputeBoundingBox(sampleRays, forcedViewSpace);

	// initialise termination criteria
	mVspTree.mTermMinProbability *= mBoundingBox.GetVolume();
	mVspTree.mGlobalCostMisses = 0;

	// get clipped rays
	mVspTree.PreprocessRays(sampleRays, rays);

	const int pvsSize = mVspTree.EvalPvsSize(rays);
	
	Debug <<  "pvs size: " << (int)pvsSize << endl;
	Debug <<  "rays size: " << (int)rays.size() << endl;

	//-- prepare view space partition

	// add first candidate for view space partition
	VspLeaf *leaf = new VspLeaf();
	mVspTree.mRoot = leaf;

	const float prop = mVspTree.mBoundingBox.GetVolume();

	// first vsp traversal data
	VspTree::VspTraversalData vData(leaf,
									0,
									&rays,
									pvsSize,	
									prop,
									mVspTree.mBoundingBox);

	// create first view cell
	mVspTree.CreateViewCell(vData, false);
		
#if WORK_WITH_VIEWCELL_PVS
	
	// add first view cell to all the objects view cell pvs
	ObjectPvsMap::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

	// compute first split candidate
	VspTree::VspSplitCandidate *splitCandidate = 
		new VspTree::VspSplitCandidate(vData);
    mVspTree.EvalSplitCandidate(*splitCandidate);

	mVspTree.mTotalCost = (float)pvsSize;
	mVspTree.EvalSubdivisionStats(*splitCandidate);

	return splitCandidate;
}


OspTree::OspSplitCandidate * HierarchyManager::PrepareOsp(const VssRayContainer &sampleRays,
														  const ObjectContainer &objects,
														  AxisAlignedBox3 *forcedObjectSpace,
														  RayInfoContainer &rays)
{
	// store pointer to this tree
	OspTree::OspSplitCandidate::sOspTree = &mOspTree;
	mOspTree.mOspStats.nodes = 1;
	
	// compute bounding box from objects
	mOspTree.ComputeBoundingBox(objects, forcedObjectSpace);

	mOspTree.mTermMinProbability *= mBoundingBox.GetVolume();
	mGlobalCostMisses = 0;

	// create new root
	KdLeaf *kdleaf = new KdLeaf(NULL, 0);
	kdleaf->mObjects = objects;
	mOspTree.mRoot = kdleaf;
	
	// get clipped rays
	mOspTree.PreprocessRays(kdleaf, sampleRays, rays);


	// probabilty is voume of all "seen" view cells
#if 1
	const float prop = mOspTree.EvalViewCellsVolume(kdleaf, rays);
#else
	const float prop = mVspTree.GetBoundingBox().GetVolume();
#endif

	//-- add first candidate for object space partition

	// create osp traversal data
	OspTree::OspTraversalData oData(kdleaf,
									0,
									&rays,
									0,//(int)objects.size(),
									prop,
									mOspTree.mBoundingBox);

		
	// first split candidate
	OspTree::OspSplitCandidate *oSplitCandidate = 
		new OspTree::OspSplitCandidate(oData);

	mOspTree.UpdateViewCellsPvs(kdleaf, rays);

	mOspTree.EvalSplitCandidate(*oSplitCandidate);

	mOspTree.mTotalCost = (float)objects.size() * prop / 
		mVspTree.GetBoundingBox().GetVolume();

	mOspTree.EvalSubdivisionStats(*oSplitCandidate);

	return oSplitCandidate;
}


void HierarchyManager::PrepareConstruction(const VssRayContainer &sampleRays,
										   const ObjectContainer &objects,
										   AxisAlignedBox3 *forcedViewSpace,
										   RayInfoContainer &viewSpaceRays,
										   RayInfoContainer &objectSpaceRays)
{
	SplitCandidate *vsc = 
		PrepareVsp(sampleRays, forcedViewSpace, viewSpaceRays);
	mTQueue.Push(vsc);

	SplitCandidate *osc = 
		PrepareOsp(sampleRays, objects, forcedViewSpace, objectSpaceRays);

	mTQueue.Push(osc);
}


void HierarchyManager::EvalSubdivisionStats(const SplitCandidate &tData)
{
	const float costDecr = tData.GetRenderCostDecrease();

	//mTotalCost -= costDecr;
	// mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;

	AddSubdivisionStats(mOspTree.mOspStats.Leaves() + mVspTree.mVspStats.Leaves(),
						costDecr,
						mTotalCost
						);
}


void HierarchyManager::AddSubdivisionStats(const int splits,
										   const float renderCostDecr,
										   const float totalRenderCost)
{
	mSubdivisionStats 
			<< "#Splits\n" << splits << endl
			<< "#RenderCostDecrease\n" << renderCostDecr << endl 
			<< "#TotalRenderCost\n" << totalRenderCost << endl;
			//<< "#AvgRenderCost\n" << avgRenderCost << endl;
}


bool HierarchyManager::GlobalTerminationCriteriaMet(SplitCandidate *candidate) const
{
	// TODO matt: make virtual by creating superclasss for traveral data
	return candidate->GlobalTerminationCriteriaMet();
}


void HierarchyManager::Construct(const VssRayContainer &sampleRays,
								 const ObjectContainer &objects,
								 AxisAlignedBox3 *forcedViewSpace)
{
	RayInfoContainer *objectSpaceRays = new RayInfoContainer();
	RayInfoContainer *viewSpaceRays = new RayInfoContainer();

	// prepare vsp and osp trees for traversal
	PrepareConstruction(sampleRays, objects, forcedViewSpace, *viewSpaceRays, *objectSpaceRays);

	mVspTree.mVspStats.Reset();
	mVspTree.mVspStats.Start();

	cout << "Constructing view space / object space tree ... \n";
	const long startTime = GetTime();	
	
	RunConstruction(true);

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

	mVspTree.mVspStats.Stop();
}


bool HierarchyManager::SubdivideSplitCandidate(SplitCandidate *sc)
{
	const bool globalTerminationCriteriaMet = 
			GlobalTerminationCriteriaMet(sc);

	const bool vspSplit = (sc->Type() == SplitCandidate::VIEW_SPACE);

	if (vspSplit)
	{
		VspNode *n = mVspTree.Subdivide(mTQueue, sc, globalTerminationCriteriaMet);
	}
	else
	{
		KdNode *n = mOspTree.Subdivide(mTQueue, sc, globalTerminationCriteriaMet);
	}
	
	return !globalTerminationCriteriaMet;
}


void HierarchyManager::RunConstruction(const bool repair)
{
	int numNodes = 0;

	while (!FinishedConstruction())
	{
		SplitCandidate *splitCandidate = NextSplitCandidate();
	    
		mTotalCost -= splitCandidate->GetRenderCostDecrease();

		// cost ratio of cost decrease / totalCost
		const float costRatio = splitCandidate->GetRenderCostDecrease() / mTotalCost;

		if (costRatio < mTermMinGlobalCostRatio)
			++ mGlobalCostMisses;
	
		/*Debug << "\n**********" << endl
			  << "total cost: " << mTotalCost << " render cost decr: " 
			  << splitCandidate->GetRenderCostDecrease() 
			  << " cost ratio: " << costRatio << endl << endl;*/

		//-- subdivide leaf node

		if (SubdivideSplitCandidate(splitCandidate))
		{
			// subdivision successful
			EvalSubdivisionStats(*splitCandidate);
		
			// reevaluate candidates affected by the split
			// for view space splits, this would be object space splits
			// and other way round
			if (repair)
				RepairQueue();

			cout << "candidate: " << splitCandidate->Type() << ", priority: " << splitCandidate->GetPriority() << endl;
		}

		DEL_PTR(splitCandidate);
	}
}


void HierarchyManager::Construct2(const VssRayContainer &sampleRays,
								  const ObjectContainer &objects,
								  AxisAlignedBox3 *forcedViewSpace)
{
	// rays clipped in view space and in object space
	RayInfoContainer *viewSpaceRays = new RayInfoContainer();
	RayInfoContainer *objectSpaceRays = new RayInfoContainer();


	/////////////////////////////////////////////////////////////
	// view space space partition
	/////////////////////////////////////////////////////////////

#if 0
	// makes no sense otherwise because only one kd cell available
	// during view space partition
	const bool savedCountMethod = mVspTree.mUseKdPvsForHeuristics;
	const bool savedStoreMethod = mVspTree.mStoreKdPvs;
	
	mVspTree.mUseKdPvsForHeuristics = false;
	mVspTree.mStoreKdPvs = false;
#endif

	VspTree::VspSplitCandidate *vsc = 
		PrepareVsp(sampleRays, forcedViewSpace, *viewSpaceRays);

	// add to queue
	mTQueue.Push(vsc);

	long startTime = GetTime();
	cout << "starting vsp contruction ... " << endl;

	mVspTree.mVspStats.Reset();
	mVspTree.mVspStats.Start();

	int i = 0;

	// all objects can be seen from everywhere
	mTotalCost = (float)vsc->mParentData.mPvs;

	const bool repairQueue = false;

	// process view space candidates
	RunConstruction(repairQueue);

	cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl;
	mVspTree.mVspStats.Stop();
	


	/////////////////////////////////////////////////////////////
	// object space partition
	/////////////////////////////////////////////////////////////


	Debug << "\n$$$$$$$$$ osp tree construction $$$$$$$$$$\n" << endl;
	cout << "starting osp contruction ... " << endl;

	// start with one big kd cell - all objects can be seen from everywhere
	// note: only true for view space = object space

	// compute first candidate
	OspTree::OspSplitCandidate *osc =
		PrepareOsp(sampleRays, objects, forcedViewSpace, *objectSpaceRays);

	Debug << "reseting cost, new total cost: " << mTotalCost << endl;
	mTotalCost = mOspTree.mTotalCost;

    mTQueue.Push(osc);

	mOspTree.mOspStats.Reset();
	mOspTree.mOspStats.Start();

	startTime = GetTime();
	
	// process object space candidates
	RunConstruction(repairQueue);
	
	cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;

	mOspTree.mOspStats.Stop();

	float rc = mOspTree.EvalRenderCost(sampleRays);

	Debug << "here47 My render cost evalulation: " << rc << endl;

#if 0
	// reset parameters
	mVspTree.mUseKdPvsForHeuristics = savedCountMethod;
	mVspTree.mStoreKdPvs = savedStoreMethod;
#endif
}


void HierarchyManager::Construct3(const VssRayContainer &sampleRays,
								  const ObjectContainer &objects,
								  AxisAlignedBox3 *forcedViewSpace)
{
	// only view space partition
	// object kd tree is taken for osp

	mVspTree.mVspStats.Reset();
	mVspTree.mVspStats.Start();

	RayInfoContainer *viewSpaceRays = new RayInfoContainer();
	
	SplitCandidate *sc = PrepareVsp(sampleRays, forcedViewSpace, *viewSpaceRays);
	mTQueue.Push(sc);

	cout << "starting vsp contruction ... " << endl;

	long startTime = GetTime();

	RunConstruction(false);

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


bool HierarchyManager::FinishedConstruction() const
{
	return mTQueue.Empty();
}


void HierarchyManager::CollectDirtyCandidates(vector<SplitCandidate *> &dirtyList)
{	
	// we have either a object space or view space split
	if (mCurrentCandidate->Type() == SplitCandidate::VIEW_SPACE)
	{
		VspTree::VspSplitCandidate *sc = 
			dynamic_cast<VspTree::VspSplitCandidate *>(mCurrentCandidate);

		mVspTree.CollectDirtyCandidates(sc, dirtyList);
	}
	else // object space split
	{			
		OspTree::OspSplitCandidate *sc = 
			dynamic_cast<OspTree::OspSplitCandidate *>(mCurrentCandidate);

		mOspTree.CollectDirtyCandidates(sc, dirtyList);
	}
}


void HierarchyManager::RepairQueue()
{
	// TODO
	// for each update of the view space partition:
	// the candidates from object space partition which
	// have been afected by the view space split (the kd split candidates
	// which saw the view cell which was split) must be reevaluated 
	// (maybe not locally, just reinsert them into the queue)
	//
	// vice versa for the view cells
	// for each update of the object space partition
	// reevaluate split candidate for view cells which saw the split kd cell
	// 
	// the priority queue update can be solved by implementing a binary heap
	// (explicit data structure, binary tree)
	// *) inserting and removal is efficient
	// *) search is not efficient => store queue position with each
	// split candidate

	// collect list of "dirty" candidates
	vector<SplitCandidate *> dirtyList;
	CollectDirtyCandidates(dirtyList);
	
	//-- reevaluate the dirty list
	vector<SubdivisionCandidate *>::const_iterator sit, sit_end = dirtyList.end();

	for (sit = dirtyList.begin(); sit != sit_end; ++ sit)
	{
		SplitCandidate* sc = *sit;

		mTQueue.Erase(sc);

		// reevaluate
		sc->EvalPriority();

		// reinsert
		mTQueue.Push(sc);
	}
}

}