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

namespace GtpVisibilityPreprocessor {

#define USE_FIXEDPOINT_T 0


//-- static members

int VspOspTree::sFrontId = 0;
int VspOspTree::sBackId = 0;
int VspOspTree::sFrontAndBackId = 0;



// 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;
	if (pvs > upper)
		return (float)upper;

	return (float)pvs;
}


int VspNode::sMailId = 1;


/******************************************************************/
/*                  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 *b, VspNode *f) 
{
    mBack = b;
    mFront = f;
}


AxisAlignedBox3 VspInterior::GetBox() const
{
	return mBox;
}


void VspInterior::SetBox(const AxisAlignedBox3 &box)
{
	mBox = box;
}


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

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


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


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

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 VspOspTree implementation                      */
/******************************************************************************/


VspOspTree::VspOspTree():
mRoot(NULL),
mOutOfBoundsCell(NULL),
mStoreRays(false),
mUseRandomAxis(false),
mTimeStamp(1)
{
	bool randomize = false;
	Environment::GetSingleton()->GetBoolValue("VspOspTree.Construction.randomize", randomize);
	if (randomize)
		Randomize(); // initialise random generator for heuristics

	//-- termination criteria for autopartition
	Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.minPvs", mTermMinPvs);
	Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.minRays", mTermMinRays);
	Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.minProbability", mTermMinProbability);
	Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxRayContribution", mTermMaxRayContribution);
	Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxCostRatio", mTermMaxCostRatio);
	Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.missTolerance", mTermMissTolerance);
	Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.maxViewCells", mMaxViewCells);

	//-- max cost ratio for early tree termination
	Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxCostRatio", mTermMaxCostRatio);

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

	// HACK//mTermMinPolygons = 25;

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

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

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

	Environment::GetSingleton()->GetBoolValue("VspOspTree.useCostHeuristics", mUseCostHeuristics);
	Environment::GetSingleton()->GetBoolValue("VspOspTree.simulateOctree", mCirculatingAxis);
	Environment::GetSingleton()->GetBoolValue("VspOspTree.useRandomAxis", mUseRandomAxis);
	Environment::GetSingleton()->GetIntValue("VspOspTree.nodePriorityQueueType", mNodePriorityQueueType);
	
	char subdivisionStatsLog[100];
	Environment::GetSingleton()->GetStringValue("VspOspTree.subdivisionStats", subdivisionStatsLog);
	mSubdivisionStats.open(subdivisionStatsLog);

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

	//-- debug output

	Debug << "******* VSP BSP 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 << "use random axis: " << mUseRandomAxis << endl;
	Debug << "priority queue type: " << mNodePriorityQueueType << endl;
	Debug << "circulating axis: " << mCirculatingAxis << endl;
	Debug << "minband: " << mMinBand << endl;
	Debug << "maxband: " << mMaxBand << endl;
	

	mSplitCandidates = new vector<SortableEntry>;

	Debug << endl;
}


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


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

	return mOutOfBoundsCell;
}


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


VspOspTree::~VspOspTree()
{
	DEL_PTR(mRoot);
	DEL_PTR(mSplitCandidates);
}


void VspOspTree::Construct(const VssRayContainer &sampleRays,
						   AxisAlignedBox3 *forcedBoundingBox)
{
	mVspStats.nodes = 1;
	mBox.Initialize(); 	// initialise BSP tree bounding box

	if (forcedBoundingBox)
		mBox = *forcedBoundingBox;
	
	PolygonContainer polys;
	RayInfoContainer *rays = new RayInfoContainer();

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

	long startTime = GetTime();

	cout << "Extracting polygons from rays ... ";

	Intersectable::NewMail();

	int numObj = 0;

	//-- store rays
	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 (mBox.GetRaySegment(hray, minT, maxT))
		{
			float len = ray->Length();

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

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

	mTermMinProbability *= mBox.GetVolume();

	mGlobalCostMisses = 0;

	cout << "finished" << endl;

	// use split cost priority queue
	Construct(rays);
}


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


void VspOspTree::Construct(RayInfoContainer *rays)
{
	VspOspSplitQueue tQueue;

	mRoot = new VspLeaf();

	const float prop = mBox.GetVolume();

	VspOspTraversalData tData(mRoot,
							  0,
							  rays,
                              ComputePvsSize(*rays),
							  prop,
							  mBox);


	// compute first split candidate
	VspOspSplitCandidate splitCandidate;
	EvalSplitCandidate(tData, splitCandidate);

	tQueue.push(splitCandidate);

	mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
	mTotalPvsSize = tData.mPvs;
	
	mSubdivisionStats 
			<< "#ViewCells\n1\n" <<  endl
			<< "#RenderCostDecrease\n0\n" << endl 
			<< "#SplitCandidateCost\n0\n" << endl
			<< "#TotalRenderCost\n" << mTotalCost << endl
			<< "#AvgRenderCost\n" << mTotalPvsSize << endl;

	Debug << "total cost: " << mTotalCost << endl;
	
	
	mVspStats.Start();
	cout << "Constructing vsp bsp tree ... \n";

	long startTime = GetTime();	
	int nLeaves = 500;
	int nViewCells = 500;

	// used for intermediate time measurements and progress
	long interTime = GetTime();	

	mOutOfMemory = false;

	mCreatedViewCells = 0;
	
	while (!tQueue.empty())
	{
		splitCandidate = tQueue.top();
	    tQueue.pop();		

		// cost ratio of cost decrease / totalCost
		float costRatio = splitCandidate.GetCost() / mTotalCost;

		//Debug << "cost ratio: " << costRatio << endl;

		if (costRatio < mTermMinGlobalCostRatio)
			++ mGlobalCostMisses;
		
		if (0 && !mOutOfMemory)
		{
			float mem = GetMemUsage();

			if (mem > mMaxMemory)
			{
				mOutOfMemory = true;
				Debug << "memory limit reached: " << mem << endl;
			}
		}

		// subdivide leaf node
		VspNode *r = Subdivide(tQueue, splitCandidate);

		if (r == mRoot)
			Debug << "VSP BSP tree construction time spent at root: "
				  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;

		if (mVspStats.Leaves() >= nLeaves)
		{
			nLeaves += 500;

			cout << "leaves=" << mVspStats.Leaves() << endl;
			Debug << "needed "
				  << TimeDiff(interTime, GetTime())*1e-3 
				  << " secs to create 500 view cells" << endl;
			interTime = GetTime();
		}

		if (mCreatedViewCells == nViewCells)
		{
			nViewCells += 500;

			cout << "generated " << mCreatedViewCells << " viewcells" << endl;
		}
	}

	Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
	cout << "finished\n";

	mVspStats.Stop();
}


bool VspOspTree::LocalTerminationCriteriaMet(const VspOspTraversalData &data) const
{
	return
		(((int)data.mRays->size() <= mTermMinRays) ||
		 (data.mPvs <= mTermMinPvs)   ||
		 (data.mProbability <= mTermMinProbability) ||
		 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
		 (data.mDepth >= mTermMaxDepth));
}


bool VspOspTree::GlobalTerminationCriteriaMet(const VspOspTraversalData &data) const
{
	return
		(mOutOfMemory || 
		(mVspStats.Leaves() >= mMaxViewCells) || 
		(mGlobalCostMisses >= mTermGlobalCostMissTolerance));
}


// subdivide using a split plane queue
VspNode *VspOspTree::Subdivide(VspOspSplitQueue &tQueue,
							   VspOspSplitCandidate &splitCandidate)
{
	VspOspTraversalData &tData = splitCandidate.mParentData;

	VspNode *newNode = tData.mNode;

	if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
	{	
		VspOspTraversalData tFrontData;
		VspOspTraversalData tBackData;

		//-- continue subdivision
		
		// create new interior node and two leaf node
		const AxisAlignedPlane splitPlane = splitCandidate.mSplitPlane;

		newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData);
	
		const int maxCostMisses = splitCandidate.mMaxCostMisses;


		// how often was max cost ratio missed in this branch?
		tFrontData.mMaxCostMisses = maxCostMisses;
		tBackData.mMaxCostMisses = maxCostMisses;
			
		
		if (1)
		{
			float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
			float cBack = (float)tBackData.mPvs * tBackData.mProbability;
			float cData = (float)tData.mPvs * tData.mProbability;

			
			float costDecr = 
				(cFront + cBack - cData) / mBox.GetVolume();

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

			mSubdivisionStats 
					<< "#ViewCells\n" << mVspStats.Leaves() << endl
					<< "#RenderCostDecrease\n" << -costDecr << endl
					<< "#SplitCandidateCost\n" << splitCandidate.GetCost() << endl
					<< "#TotalRenderCost\n" << mTotalCost << endl
					<< "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mVspStats.Leaves() << endl;
		}

	
		//-- push the new split candidates on the stack
		VspOspSplitCandidate frontCandidate;
		VspOspSplitCandidate backCandidate;

		EvalSplitCandidate(tFrontData, frontCandidate);
		EvalSplitCandidate(tBackData, backCandidate);
	
		tQueue.push(frontCandidate);
		tQueue.push(backCandidate);
	
		// delete old leaf node
		DEL_PTR(tData.mNode);
	}


	//-- terminate traversal and create new view cell
	if (newNode->IsLeaf())
	{
		VspLeaf *leaf = dynamic_cast<VspLeaf *>(newNode);
		VspViewCell *viewCell = new VspViewCell();

        leaf->SetViewCell(viewCell);
		
		//-- update pvs
		int conSamp = 0;
		float sampCon = 0.0f;
		AddToPvs(leaf, *tData.mRays, sampCon, conSamp);

		// update scalar pvs size value
		mViewCellsManager->SetScalarPvsSize(viewCell, viewCell->GetPvs().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);
			}
		}

		// should I check here?
		viewCell->mLeaf = leaf;

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

		EvaluateLeafStats(tData);		
	}

	//-- cleanup
	tData.Clear();

	return newNode;
}


void VspOspTree::EvalSplitCandidate(VspOspTraversalData &tData,
									VspOspSplitCandidate &splitData)
{
	VspOspTraversalData frontData;
	VspOspTraversalData backData;

	VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);

	// compute locally best split plane
	const bool success = 
		SelectPlane(tData, splitData.mSplitPlane,
					frontData.mProbability, backData.mProbability);

	//TODO
	// compute global decrease in render cost
	splitData.mRenderCost = EvalRenderCostDecrease(splitData.mSplitPlane, tData);
	splitData.mParentData = tData;
	splitData.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
}


VspInterior *VspOspTree::SubdivideNode(const AxisAlignedPlane &splitPlane,
									   VspOspTraversalData &tData,
									   VspOspTraversalData &frontData,
									   VspOspTraversalData &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;
	frontData.mRays = new RayInfoContainer();
	
	backData.mDepth = tData.mDepth + 1;
	backData.mRays = new RayInfoContainer();
	
	//-- subdivide rays
	SplitRays(splitPlane,
			  *tData.mRays,
			  *frontData.mRays,
			  *backData.mRays);

	//-- compute pvs
	frontData.mPvs = ComputePvsSize(*frontData.mRays);
	backData.mPvs = ComputePvsSize(*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;
	}

	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);
	}
	else // new root
	{
		mRoot = interior;
	}

	// and setup child links
	interior->SetupChildLinks(new VspLeaf(interior), new VspLeaf(interior));

	frontData.mNode = interior->GetFront();
	backData.mNode = interior->GetBack();

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


void VspOspTree::AddToPvs(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;

		if (ray->mTerminationObject) 
		{
			if (vc->AddPvsSample(ray->mTerminationObject, ray->mPdf, contribution))
				madeContrib = true;
			sc += contribution;
		}
	  
		if (ray->mOriginObject) 
		{
			if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
				madeContrib = true;
			sc += contribution;
		}
	  
	  sampleContributions += sc;

	  if (madeContrib)
		  ++ contributingSamples;
		
	  //leaf->mVssRays.push_back(ray);
	}
}


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

	int requestedSize = 2 * (int)(rays.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;

	// float values => don't compare with exact values
	if (0)
	{
		minBand += Limits::Small;
		maxBand -= Limits::Small;
	}

	// insert all queries
	for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
	{
		const bool positive = (*ri).mRay->HasPosDir(axis);
				
		pos = (*ri).ExtrapOrigin(axis);
		// clamp to min / max band
		if (0) ClipValue(pos, minBand, maxBand);
		
		mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax, 
									pos, (*ri).mRay));

		pos = (*ri).ExtrapTermination(axis);
		// clamp to min / max band
		if (0) ClipValue(pos, minBand, maxBand);

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

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


float VspOspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
										  const AxisAlignedBox3 &box,
										  const int pvsSize,
										  const int axis,
                                          float &position)
{
	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(rays, axis, minBand, maxBand);

	// 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 border can be found, take mid split
	position = minBox + 0.5f * sizeBox;
	
	// the relative cost ratio
	float ratio = /*Limits::Infinity;*/99999999.0f;
	bool splitPlaneFound = false;

	Intersectable::NewMail();

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

	// set all object as belonging to the front pvs
	for(ri = rays.begin(); ri != ri_end; ++ ri)
	{
		Intersectable *oObject = (*ri).mRay->mOriginObject;
		Intersectable *tObject = (*ri).mRay->mTerminationObject;

		if (oObject)
		{
			if (!oObject->Mailed())
			{
				oObject->Mail();
				oObject->mCounter = 1;
			}
			else
			{
				++ oObject->mCounter;
			}
		}
		if (tObject)
		{
			if (!tObject->Mailed())
			{
				tObject->Mail();
				tObject->mCounter = 1;
			}
			else
			{
				++ tObject->mCounter;
			}
		}
	}

	Intersectable::NewMail();

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

	for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
	{
		VssRay *ray;
		ray = (*ci).ray;
		
		Intersectable *oObject = ray->mOriginObject;
		Intersectable *tObject = ray->mTerminationObject;
		

		switch ((*ci).type)
		{
			case SortableEntry::ERayMin:
				{
					if (oObject && !oObject->Mailed())
					{
						oObject->Mail();
						++ pvsl;
					}
					if (tObject && !tObject->Mailed())
					{
						tObject->Mail();
						++ pvsl;
					}
					break;
				}
			case SortableEntry::ERayMax:
				{
					if (oObject)
					{
						if (-- oObject->mCounter == 0)
							-- pvsr;
					}

					if (tObject)
					{
						if (-- tObject->mCounter == 0)
							-- pvsr;
					}

					break;
				}
		}

		
		
		// Note: sufficient to compare size of bounding boxes of front and back side?
		if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
		{
			sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);

			//Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << endl;
			//Debug << "cost= " << sum << endl;

			if (sum < minSum)
			{
				splitPlaneFound = true;

				minSum = sum;
				position = (*ci).value;
				
				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;
	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

	if (splitPlaneFound)
	{
		ratio = newRenderCost / (oldRenderCost + Limits::Small);
	}
	//if (axis != 1)
	//Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
	 //    <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;

	return ratio;
}


float VspOspTree::SelectPlane(const VspOspTraversalData &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 || mUseRandomAxis || mCirculatingAxis;

	if (mUseRandomAxis)
		sAxis = Random(3);
	else if (mCirculatingAxis)
		sAxis = (tData.mAxis + 1) % 3;
	else if (mOnlyDrivingAxis)
		sAxis = box.Size().DrivingAxis();


	for (int axis = 0; axis < 3 ; ++ axis)
	{
		if (!useSpecialAxis || (axis == sAxis))
		{
			//-- place split plane using heuristics

			if (mUseCostHeuristics)
			{
				nCostRatio[axis] =
					BestCostRatioHeuristics(*tData.mRays,
										    box,
											tData.mPvs,
											axis,
											nPosition[axis]);			
			}
			else //-- split plane position is spatial median
			{
				nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;

				nCostRatio[axis] = EvalSplitCost(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;
	pFront = nProbFront[bestAxis];
	pBack = nProbBack[bestAxis];

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

	return nCostRatio[bestAxis];
}


inline void VspOspTree::GenerateUniqueIdsForPvs()
{
	Intersectable::NewMail(); sBackId = Intersectable::sMailId;
	Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
	Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
}


float VspOspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
										 const VspOspTraversalData &data) 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;


	// create unique ids for pvs heuristics
	GenerateUniqueIdsForPvs();
	
	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);

		// find front and back pvs for origing and termination object
		AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
		AddObjToPvs(ray->mOriginObject, 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(totalPvs, lowerPvsLimit, upperPvsLimit);
    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
			
	const float oldRenderCost = pOverall * penaltyOld;
	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

	//Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
	const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBox.GetVolume();
	

#if 1
	// take render cost of node into account to avoid being stuck in a local minimum
	const float normalizedOldRenderCost = oldRenderCost / mBox.GetVolume();
	return 0.99f * renderCostDecrease + 0.01f * normalizedOldRenderCost;
#else
	return renderCostDecrease;
#endif
}



float VspOspTree::EvalSplitCost(const VspOspTraversalData &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
	GenerateUniqueIdsForPvs();

	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);
	
		AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
		AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
	}

	//-- pvs heuristics
	float pOverall;

	//-- compute heurstics
	//-- we take simplified computation for mid split
		
	pOverall = data.mProbability;
	pBack = pFront = pOverall * 0.5f;
	
	
#ifdef _DEBUG
	Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
	Debug << pFront << " " << pBack << " " << pOverall << endl;
#endif

	
	const float newCost = pvsBack * pBack + pvsFront * pFront;
	const float oldCost = (float)pvsSize * pOverall + Limits::Small;

	return  (mCtDivCi + newCost) / oldCost;
}


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

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

	// new object
	if ((obj->mMailbox != sFrontId) &&
		(obj->mMailbox != sBackId) &&
		(obj->mMailbox != sFrontAndBackId))
	{
		totalPvs += renderCost;
	}

	// TODO: does this really belong to no pvs?
	//if (cf == Ray::COINCIDENT) return;

	// object belongs to both PVS
	if (cf >= 0)
	{
		if ((obj->mMailbox != sFrontId) &&
			(obj->mMailbox != sFrontAndBackId))
		{
			frontPvs += renderCost;
		
			if (obj->mMailbox == sBackId)
				obj->mMailbox = sFrontAndBackId;
			else
				obj->mMailbox = sFrontId;
		}
	}

	if (cf <= 0)
	{
		if ((obj->mMailbox != sBackId) &&
			(obj->mMailbox != sFrontAndBackId))
		{
			backPvs += renderCost;
		
			if (obj->mMailbox == sFrontId)
				obj->mMailbox = sFrontAndBackId;
			else
				obj->mMailbox = sBackId;
		}
	}
}


void VspOspTree::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().GetSize() <= maxPvsSize)))
			{
				leaves.push_back(leaf);
			}
		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);

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


AxisAlignedBox3 VspOspTree::GetBoundingBox() const
{
	return mBox;
}


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


void VspOspTree::EvaluateLeafStats(const VspOspTraversalData &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 << "), "
	  //		  << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
		  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
		  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
		  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
#endif
}


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


void VspOspTree::CollapseViewCells()
{
// TODO
#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 VspOspTree::CollectRays(VssRayContainer &rays)
{
	vector<VspLeaf *> 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 VspOspTree::ValidateTree()
{
	stack<VspNode *> nodeStack;

	if (!mRoot)
		return;

	nodeStack.push(mRoot);
	
	mVspStats.invalidLeaves = 0;
	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 VspOspTree::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 VspOspTree::FindNeighbors(VspLeaf *n,
							  vector<VspLeaf *> &neighbors,
							  const bool onlyUnmailed) const
{
	stack<VspNode *> nodeStack;

	nodeStack.push(mRoot);

	const AxisAlignedBox3 box = GetBBox(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 *VspOspTree::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 (GetBBox(interior->GetBack()).Side(plane) < 0)
				next = interior->GetFront();
            else
				if (GetBBox(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 *VspOspTree::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;
}


int VspOspTree::ComputePvsSize(const RayInfoContainer &rays) const
{
	int pvsSize = 0;

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

	Intersectable::NewMail();

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

	return pvsSize;
}


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


int VspOspTree::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();

			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 VspOspTree::CastRay(Ray &ray)
{
	int hits = 0;

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

	float maxt, mint;

	if (!mBox.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;
}


VspNode *VspOspTree::CollapseTree(VspNode *node, int &collapsed)
{
// TODO
#if HAS_TO_BE_REDONE
	if (node->IsLeaf())
		return node;

	VspInterior *interior = dynamic_cast<VspInterior *>(node);

	VspNode *front = CollapseTree(interior->GetFront(), collapsed);
	VspNode *back = CollapseTree(interior->GetBack(), collapsed);

	if (front->IsLeaf() && back->IsLeaf())
	{
		VspLeaf *frontLeaf = dynamic_cast<VspLeaf *>(front);
		VspLeaf *backLeaf = dynamic_cast<VspLeaf *>(back);

		//-- collapse tree
		if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
		{
			BspViewCell *vc = frontLeaf->GetViewCell();

			VspLeaf *leaf = new VspLeaf(interior->GetParent(), vc);
			leaf->SetTreeValid(frontLeaf->TreeValid());

			// replace a link from node's parent
			if (leaf->GetParent())
				leaf->GetParent()->ReplaceChildLink(node, leaf);
			else
				mRoot = leaf;

			++ collapsed;
			delete interior;

			return leaf;
		}
	}
#endif
	return node;
}


int VspOspTree::CollapseTree()
{
	int collapsed = 0;
	//TODO
#if HAS_TO_BE_REDONE
	(void)CollapseTree(mRoot, collapsed);

	// revalidate leaves
	RepairViewCellsLeafLists();
#endif
	return collapsed;
}


void VspOspTree::RepairViewCellsLeafLists()
{
// TODO
#if HAS_TO_BE_REDONE
	// list not valid anymore => clear
	stack<VspNode *> nodeStack;
	nodeStack.push(mRoot);

	ViewCell::NewMail();

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

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

			BspViewCell *viewCell = leaf->GetViewCell();

			if (!viewCell->Mailed())
			{
				viewCell->mLeaves.clear();
				viewCell->Mail();
			}
	
			viewCell->mLeaves.push_back(leaf);

		}
		else
		{
			VspInterior *interior = dynamic_cast<VspInterior *>(node);

			nodeStack.push(interior->GetFront());
			nodeStack.push(interior->GetBack());
		}
	}
// TODO
#endif
}



ViewCell *VspOspTree::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 VspOspTree::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 VspOspTree::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());
		}
	}
}

#if ZIPPED_VIEWCELLS
bool VspOspTree::Export(ogzstream &stream)
#else
bool VspOspTree::Export(ofstream &stream)
#endif
{
	ExportNode(mRoot, stream);

	return true;
}


#if ZIPPED_VIEWCELLS
void VspOspTree::ExportNode(VspNode *node, ogzstream &stream)
#else
void VspOspTree::ExportNode(VspNode *node, ofstream &stream)
#endif
{
	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 VspOspTree::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 = plane.ComputeRayIntersection(bRay, 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 VspOspTree::GetBBox(VspNode *node) const
{
	if (!node->GetParent())
		return mBox;

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

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

	AxisAlignedBox3 box(parent->GetBox());

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

	return box;
}


}