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

#include "Plane3.h"
#include "VspBspTree.h"
#include "Mesh.h"
#include "common.h"
#include "ViewCell.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include "Exporter.h"
#include "Plane3.h"
#include "ViewCellBsp.h"
#include "ViewCellsManager.h"
#include "Beam.h"

#define USE_FIXEDPOINT_T 0


//-- static members

/** Evaluates split plane classification with respect to the plane's
	contribution for a minimum number of ray splits.
*/
const float VspBspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
/** Evaluates split plane classification with respect to the plane's
	contribution for  balanced rays.
*/
const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};


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




// pvs penalty can be different from pvs size
inline 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;
}




/******************************************************************************/
/*                       class VspBspTree implementation                      */
/******************************************************************************/


VspBspTree::VspBspTree():
mRoot(NULL),
mUseAreaForPvs(false),
mCostNormalizer(Limits::Small),
mViewCellsManager(NULL),
mOutOfBoundsCell(NULL),
mStoreRays(false),
mRenderCostWeight(0.5),
mUseRandomAxis(false),
mTimeStamp(0)
{
	bool randomize = false;
	environment->GetBoolValue("VspBspTree.Construction.randomize", randomize);
	if (randomize)
		Randomize(); // initialise random generator for heuristics

	//-- termination criteria for autopartition
	environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
	environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
	environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
	environment->GetFloatValue("VspBspTree.Termination.minProbability", mTermMinProbability);
	environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
	environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
	environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
	environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
	environment->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
	//-- max cost ratio for early tree termination
	environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);

	// HACK
	mTermMinPolygons = 25;

	//-- factors for bsp tree split plane heuristics
	environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
	environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
	environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);


	//-- partition criteria
	environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
	environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
	environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);

	environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
	environment->GetIntValue("VspBspTree.maxTests", mMaxTests);

	// if only the driving axis is used for axis aligned split
	environment->GetBoolValue("VspBspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
	
	//-- termination criteria for axis aligned split
	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxRayContribution", 
		mTermMaxRayContriForAxisAligned);
	environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
		mTermMinRaysForAxisAligned);

	//environment->GetFloatValue("VspBspTree.maxTotalMemory", mMaxTotalMemory);
	environment->GetFloatValue("VspBspTree.maxStaticMemory", mMaxMemory);

	environment->GetFloatValue("VspBspTree.Construction.renderCostWeight", mRenderCostWeight);

	environment->GetBoolValue("VspBspTree.usePolygonSplitIfAvailable", mUsePolygonSplitIfAvailable);

	mSubdivisionStats.open("subdivisionStats.log");

	environment->GetBoolValue("VspBspTree.useCostHeuristics", mUseCostHeuristics);


	//-- 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 << "VSP BSP mininam accumulated ray lenght: ", mTermMinAccRayLength) << endl;
	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	Debug << "miss tolerance: " << mTermMissTolerance << endl;
	Debug << "max view cells: " << mMaxViewCells << endl;
	Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
	Debug << "max plane candidates: " << mMaxRayCandidates << endl;
	Debug << "randomize: " << randomize << endl;

	Debug << "using area for pvs: " << mUseAreaForPvs << endl;
	Debug << "render cost weight: " << mRenderCostWeight << endl;
	Debug << "Split plane strategy: ";

	if (mSplitPlaneStrategy & RANDOM_POLYGON)
	{
		Debug << "random polygon ";
	}
	if (mSplitPlaneStrategy & AXIS_ALIGNED)
	{
		Debug << "axis aligned ";
	}
	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
	{
		mCostNormalizer += mLeastRaySplitsFactor;
		Debug << "least ray splits ";
	}
	if (mSplitPlaneStrategy & BALANCED_RAYS)
	{
		mCostNormalizer += mBalancedRaysFactor;
		Debug << "balanced rays ";
	}
	if (mSplitPlaneStrategy & PVS)
	{
		mCostNormalizer += mPvsFactor;
		Debug << "pvs";
	}


	mSplitCandidates = new vector<SortableEntry>;

	Debug << endl;
}


BspViewCell *VspBspTree::GetOutOfBoundsCell()
{
	return mOutOfBoundsCell;
}


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

	return mOutOfBoundsCell;
}


const BspTreeStatistics &VspBspTree::GetStatistics() const
{
	return mBspStats;
}


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


int VspBspTree::AddMeshToPolygons(Mesh *mesh,
								  PolygonContainer &polys,
								  MeshInstance *parent)
{
	FaceContainer::const_iterator fi;

	// copy the face data to polygons
	for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
	{
		Polygon3 *poly = new Polygon3((*fi), mesh);

		if (poly->Valid(mEpsilon))
		{
			poly->mParent = parent; // set parent intersectable
			polys.push_back(poly);
		}
		else
			DEL_PTR(poly);
	}
	return (int)mesh->mFaces.size();
}


int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells,
							  PolygonContainer &polys,
							  int maxObjects)
{
	int limit = (maxObjects > 0) ?
		Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();

	int polysSize = 0;

	for (int i = 0; i < limit; ++ i)
	{
		if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
		{
			mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
			polysSize +=
				AddMeshToPolygons(viewCells[i]->GetMesh(),
								  polys,
								  viewCells[i]);
		}
	}

	return polysSize;
}


int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects,
								 PolygonContainer &polys,
								 int maxObjects)
{
	int limit = (maxObjects > 0) ?
		Min((int)objects.size(), maxObjects) : (int)objects.size();

	for (int i = 0; i < limit; ++i)
	{
		Intersectable *object = objects[i];//*it;
		Mesh *mesh = NULL;

		switch (object->Type()) // extract the meshes
		{
		case Intersectable::MESH_INSTANCE:
			mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
			break;
		case Intersectable::VIEW_CELL:
			mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
			break;
			// TODO: handle transformed mesh instances
		default:
			Debug << "intersectable type not supported" << endl;
			break;
		}

        if (mesh) // copy the mesh data to polygons
		{
			mBox.Include(object->GetBox()); // add to BSP tree aabb
			AddMeshToPolygons(mesh, polys, NULL);
		}
	}

	return (int)polys.size();
}


void VspBspTree::Construct(const VssRayContainer &sampleRays,
						   AxisAlignedBox3 *forcedBoundingBox)
{
	mBspStats.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;

	//-- extract polygons intersected by the rays
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = *rit;

		if ((mBox.IsInside(ray->mTermination) || !forcedBoundingBox) &&
			ray->mTerminationObject && 
			!ray->mTerminationObject->Mailed())
		{
			ray->mTerminationObject->Mail();
			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
			AddMeshToPolygons(obj->GetMesh(), polys, obj);
			++ numObj;
			//-- compute bounding box
			if (!forcedBoundingBox)
				mBox.Include(ray->mTermination);
		}

		if ((mBox.IsInside(ray->mOrigin) || !forcedBoundingBox) &&
			ray->mOriginObject && 
			!ray->mOriginObject->Mailed())
		{
			ray->mOriginObject->Mail();
			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
			AddMeshToPolygons(obj->GetMesh(), polys, obj);
			++ numObj;

			//-- compute bounding box
			if (!forcedBoundingBox)
				mBox.Include(ray->mOrigin);
		}
	}
	
	Debug << "maximal pvs (i.e., pvs still considered as valid) : " << mViewCellsManager->GetMaxPvsSize() << endl;

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

	// normalize
	if (mUseAreaForPvs)
		mTermMinProbability *= mBox.SurfaceArea(); 
	else 
		mTermMinProbability *= mBox.GetVolume();

	mBspStats.polys = (int)polys.size();

	cout << "finished" << endl;

	Debug << "\nPolygon extraction: " << (int)polys.size() << " polys extracted from "
		  << (int)sampleRays.size() << " rays in "
		  << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;

	Construct(polys, rays);

	// clean up polygons
	CLEAR_CONTAINER(polys);
}


// TODO: return memory usage in MB
float VspBspTree::GetMemUsage(/*const VspBspTraversalQueue &tstack*/) const
{
	return
		(sizeof(VspBspTree) +
		 (float)mBspStats.Leaves() * sizeof(BspLeaf) + 
		 // the nodes in the stack is the minimal additional number of leaves
		 //(float)tstack.size() * sizeof(BspLeaf) +
		 mBspStats.Interior() * sizeof(BspInterior) +
		 mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
}



void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
{
	VspBspTraversalQueue tQueue;

	mRoot = new BspLeaf();

	// constrruct root node geometry
	BspNodeGeometry *geom = new BspNodeGeometry();
	ConstructGeometry(mRoot, *geom);

	const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();

	VspBspTraversalData tData(mRoot,
							  new PolygonContainer(polys),
							  0,
							  rays,
                              ComputePvsSize(*rays),
							  prop,
							  geom);
#if OCTREE_HACK
tData.mAxis = 0;
#endif
	// first node is kd node, i.e. an axis aligned box
	if (1)
        tData.mIsKdNode = true;
	else
		tData.mIsKdNode = false;

	tQueue.push(tData);


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

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

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

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

	mOutOfMemory = false;

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

		if (0 && !mOutOfMemory)
		{
			float mem = GetMemUsage();

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

		// subdivide leaf node
		BspNode *r = Subdivide(tQueue, tData);

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

		if (mBspStats.Leaves() == nLeaves)
		{
			nLeaves += 500;

			cout << "leaves=" << mBspStats.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";

	mBspStats.Stop();
}


bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
{
	return
		(((int)data.mRays->size() <= mTermMinRays) ||
		 (data.mPvs <= mTermMinPvs)   ||
		 (data.mProbability <= mTermMinProbability) ||
		 (mBspStats.Leaves() >= mMaxViewCells) ||
#if 0
		 (((int)data.mPolygons->size() <= mTermMinPolygons) && !data.mPolygons->empty())||
#endif
		 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
		 (data.mDepth >= mTermMaxDepth));
}


BspNode *VspBspTree::Subdivide(VspBspTraversalQueue &tQueue,
							   VspBspTraversalData &tData)
{
	BspNode *newNode = tData.mNode;

	if (!mOutOfMemory && !TerminationCriteriaMet(tData))
	{
		PolygonContainer coincident;

		VspBspTraversalData tFrontData;
		VspBspTraversalData tBackData;

#if OCTREE_HACK
		//Debug << "new axis:" << (tData.mAxis + 1) % 3 << endl;
		tFrontData.mAxis = (tData.mAxis + 1) % 3;
		tBackData.mAxis = (tData.mAxis + 1) % 3;
#endif
		// create new interior node and two leaf nodes
		// or return leaf as it is (if maxCostRatio missed)
		newNode = SubdivideNode(tData, tFrontData, tBackData, coincident);

		if (!newNode->IsLeaf()) //-- continue subdivision
		{
			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" << mBspStats.Leaves() << endl
						<< "#RenderCostDecrease\n" << -costDecr << endl 
						<< "#TotalRenderCost\n" << mTotalCost << endl
						<< "#AvgRenderCost\n" << mTotalPvsSize / mBspStats.Leaves() << endl;
			}

			// push the children on the stack
			tQueue.push(tFrontData);
			tQueue.push(tBackData);

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

	//-- terminate traversal and create new view cell
	if (newNode->IsLeaf())
	{
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);
		BspViewCell *viewCell = new BspViewCell();
		
		leaf->SetViewCell(viewCell);
	
		//-- update pvs
		int conSamp = 0;
		float sampCon = 0.0f;
		AddToPvs(leaf, *tData.mRays, sampCon, conSamp);

		mBspStats.contributingSamples += conSamp;
		mBspStats.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?
		if (0 && !mViewCellsManager->CheckValidity(viewCell, 0, mViewCellsManager->GetMaxPvsSize()))
		{
			viewCell->SetValid(false);
			leaf->SetTreeValid(false);
			PropagateUpValidity(leaf);

			++ mBspStats.invalidLeaves;
		}
		
        viewCell->mLeaf = leaf;

		if (mUseAreaForPvs)
			viewCell->SetArea(tData.mProbability);
		else
			viewCell->SetVolume(tData.mProbability);

		leaf->mProbability = tData.mProbability;

		EvaluateLeafStats(tData);		
	}

	//-- cleanup
	tData.Clear();

	return newNode;
}


BspNode *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
								   VspBspTraversalData &frontData,
								   VspBspTraversalData &backData,
								   PolygonContainer &coincident)
{
	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
	
	// select subdivision plane
	Plane3 splitPlane;
	
	int maxCostMisses = tData.mMaxCostMisses;

	int splitAxis;

	const bool success = 
		SelectPlane(splitPlane, leaf, tData, frontData, backData, splitAxis);

	if (!success)
	{
		++ maxCostMisses;

		if (maxCostMisses > mTermMissTolerance)
		{
			// terminate branch because of max cost
			++ mBspStats.maxCostNodes;
            return leaf;
		}
	}

	//! error also computed if cost ratio is missed
	if (splitAxis < 3)
		++ mBspStats.splits[splitAxis];
	else
		++ mBspStats.polySplits;

	mBspStats.nodes += 2;

	//-- subdivide further
	BspInterior *interior = new BspInterior(splitPlane);

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

	//-- the front and back traversal data is filled with the new values
	frontData.mDepth = tData.mDepth + 1;
	frontData.mPolygons = new PolygonContainer();
	frontData.mRays = new RayInfoContainer();
	
	backData.mDepth = tData.mDepth + 1;
	backData.mPolygons = new PolygonContainer();
	backData.mRays = new RayInfoContainer();
	
	// subdivide rays
	SplitRays(interior->GetPlane(),
			  *tData.mRays,
			  *frontData.mRays,
			  *backData.mRays);

	// subdivide polygons
	SplitPolygons(interior->GetPlane(),
				  *tData.mPolygons,
	              *frontData.mPolygons,
				  *backData.mPolygons,
				  coincident);


	// how often was max cost ratio missed in this branch?
	frontData.mMaxCostMisses = maxCostMisses;
	backData.mMaxCostMisses = maxCostMisses;

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

	// split front and back node geometry and compute area
	
	// if geometry was not already computed
	if (!frontData.mGeometry && !backData.mGeometry)
	{
		frontData.mGeometry = new BspNodeGeometry();
		backData.mGeometry = new BspNodeGeometry();

		tData.mGeometry->SplitGeometry(*frontData.mGeometry,
									   *backData.mGeometry,
									   interior->GetPlane(),
									   mBox,
									   mEpsilon);
		
		if (mUseAreaForPvs)
		{
			frontData.mProbability = frontData.mGeometry->GetArea();
			backData.mProbability = backData.mGeometry->GetArea();
		}
		else
		{
			frontData.mProbability = frontData.mGeometry->GetVolume();
			backData.mProbability =  tData.mProbability - frontData.mProbability;
		}
	}
	

	//-- create front and back leaf

	BspInterior *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 BspLeaf(interior), new BspLeaf(interior));

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

	frontData.mNode->mTimeStamp = mTimeStamp;
	backData.mNode->mTimeStamp = mTimeStamp ++;

	//DEL_PTR(leaf);
	return interior;
}


void VspBspTree::AddToPvs(BspLeaf *leaf,
						  const RayInfoContainer &rays,
						  float &sampleContributions,
						  int &contributingSamples)
{
  sampleContributions = 0;
  contributingSamples = 0;
  
  RayInfoContainer::const_iterator it, it_end = rays.end();
  
  ViewCell *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->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf, contribution))
		  madeContrib = true;
		sc += contribution;
	  }
	  
	  if (ray->mOriginObject) {
		if (vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf, contribution))
		  madeContrib = true;
		sc += contribution;
	  }
	  
	  sampleContributions += sc;
	  if (madeContrib)
		  ++ contributingSamples;
		
	  //leaf->mVssRays.push_back(ray);
	}
}


void VspBspTree::SortSplitCandidates(const RayInfoContainer &rays, const int axis)
{
	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);

	// insert all queries
	for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
	{
		const bool positive = (*ri).mRay->HasPosDir(axis);
		mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
												  (*ri).ExtrapOrigin(axis), (*ri).mRay));

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

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


float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
										  const AxisAlignedBox3 &box,
										  const int pvsSize,
										  const int &axis,
                                          float &position)
{
	SortSplitCandidates(rays, axis);

	// 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, pvsr = pvsSize;

	int pvsBack = pvsl;
	int pvsFront = pvsr;

	float minBox = box.Min(axis);
	float maxBox = box.Max(axis);
	float sizeBox = maxBox - minBox;

	float minBand = minBox + 0.1f * (maxBox - minBox);
	float maxBand = minBox + 0.9f * (maxBox - minBox);

	float sum = (float)pvsSize * sizeBox;
	float minSum = 1e20f;

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

			//  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
			// cout<<"cost= "<<sum<<endl;

			if (sum < minSum)
			{
				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;

	float ratio = mPvsFactor * newRenderCost / (oldRenderCost + Limits::Small);

	//Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
	//     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;

	return ratio;
}


float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
										 const VspBspTraversalData &tData,
										 int &axis,
										 BspNodeGeometry **frontGeom,
										 BspNodeGeometry **backGeom,
										 float &pFront,
										 float &pBack,
										 const bool useKdSplit)
{
	float nPosition[3];
	float nCostRatio[3];
	float nProbFront[3];
	float nProbBack[3];

	BspNodeGeometry *nFrontGeom[3];
	BspNodeGeometry *nBackGeom[3];

	for (int i = 0; i < 3; ++ i)
	{
		nFrontGeom[i] = NULL;
		nBackGeom[i] = NULL;
	}

	int bestAxis = -1;

	// create bounding box of node geometry
	AxisAlignedBox3 box;
	box.Initialize();
	
	//TODO: for kd split geometry already is box => only take minmax vertices
	if (1)
	{
		PolygonContainer::const_iterator it, it_end = tData.mGeometry->mPolys.end();

		for(it = tData.mGeometry->mPolys.begin(); it < it_end; ++ it)
			box.Include(*(*it));
	}
	else
	{
		RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();

		for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
			box.Include((*ri).ExtrapTermination());
	}
#if OCTREE_HACK
	//Debug << "choosing axis:" << tData.mAxis << endl;
	const int sAxis = tData.mAxis;
#else
	const int sAxis = mUseRandomAxis ? Random(3) : box.Size().DrivingAxis();
#endif
	for (axis = 0; axis < 3; ++ axis)
	{
		if (!mOnlyDrivingAxis || (axis == sAxis))
		{
			if (!mUseCostHeuristics)
			{
				nFrontGeom[axis] = new BspNodeGeometry();
				nBackGeom[axis] = new BspNodeGeometry();

				nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
				Vector3 normal(0,0,0); normal[axis] = 1.0f;

				// allows faster split because we have axis aligned kd tree boxes
				if (useKdSplit)
				{
					nCostRatio[axis] = EvalAxisAlignedSplitCost(tData,
																box,
																axis,
																nPosition[axis],
																nProbFront[axis], 
																nProbBack[axis]);
					
					Vector3 pos;
					
					pos = box.Max(); pos[axis] = nPosition[axis];
					AxisAlignedBox3 bBox(box.Min(), pos);
					bBox.ExtractPolys(nBackGeom[axis]->mPolys);
					
					pos = box.Min(); pos[axis] = nPosition[axis];
					AxisAlignedBox3 fBox(pos, box.Max());
					fBox.ExtractPolys(nFrontGeom[axis]->mPolys);
				}
				else
				{
					nCostRatio[axis] = 
						EvalSplitPlaneCost(Plane3(normal, nPosition[axis]), 
										   tData, *nFrontGeom[axis], *nBackGeom[axis],
										   nProbFront[axis], nProbBack[axis]);
				}
			}
			else
			{
				nCostRatio[axis] =
					BestCostRatioHeuristics(*tData.mRays,
										    box,
											tData.mPvs,
											axis,
											nPosition[axis]);
			}

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

		}
	}

	//-- assign values
	axis = bestAxis;
	pFront = nProbFront[bestAxis];
	pBack = nProbBack[bestAxis];

	// assign best split nodes geometry 
	*frontGeom = nFrontGeom[bestAxis];
	*backGeom = nBackGeom[bestAxis];

	// and delete other geometry
	DEL_PTR(nFrontGeom[(bestAxis + 1) % 3]);
	DEL_PTR(nBackGeom[(bestAxis + 2) % 3]);

	//-- split plane
    Vector3 normal(0,0,0); normal[bestAxis] = 1;
	plane = Plane3(normal, nPosition[bestAxis]);

	return nCostRatio[bestAxis];
}


bool VspBspTree::SelectPlane(Plane3 &bestPlane,
							 BspLeaf *leaf,
							 VspBspTraversalData &data,
							 VspBspTraversalData &frontData,
							 VspBspTraversalData &backData,
							 int &splitAxis)
{
	// simplest strategy: just take next polygon
	if (mSplitPlaneStrategy & RANDOM_POLYGON)
	{
        if (!data.mPolygons->empty())
		{
			const int randIdx = 
				(int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
			Polygon3 *nextPoly = (*data.mPolygons)[randIdx];

			bestPlane = nextPoly->GetSupportingPlane();
			return true;
		}
	}

	//-- use heuristics to find appropriate plane

	// intermediate plane
	Plane3 plane;
	float lowestCost = MAX_FLOAT;
	
	// decides if the first few splits should be only axisAligned
	const bool onlyAxisAligned  = 
		(mSplitPlaneStrategy & AXIS_ALIGNED) &&
		((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
		((int)data.GetAvgRayContribution() < mTermMaxRayContriForAxisAligned);
	
	const int limit = onlyAxisAligned ? 0 : 
		Min((int)data.mPolygons->size(), mMaxPolyCandidates);

	float candidateCost;

	int maxIdx = (int)data.mPolygons->size();

	for (int i = 0; i < limit; ++ i)
	{
		// the already taken candidates are stored behind maxIdx
		// => assure that no index is taken twice
		const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
		Polygon3 *poly = (*data.mPolygons)[candidateIdx];

		// swap candidate to the end to avoid testing same plane
		std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
		//Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];

		// evaluate current candidate
		BspNodeGeometry fGeom, bGeom;
		float fArea, bArea;
		plane = poly->GetSupportingPlane();
		candidateCost = EvalSplitPlaneCost(plane, data, fGeom, bGeom, fArea, bArea);
		
		if (candidateCost < lowestCost)
		{
			bestPlane = plane;
			lowestCost = candidateCost;
		}
	}

	// cost ratio miss
	if (mUsePolygonSplitIfAvailable && !data.mPolygons->empty())
	{
		frontData.mIsKdNode = backData.mIsKdNode = false;
		if (lowestCost > mTermMaxCostRatio)
			return false;

		return true;
	}

	//-- evaluate axis aligned splits
	int axis;
	BspNodeGeometry *fGeom, *bGeom;
	float pFront, pBack;

	candidateCost = SelectAxisAlignedPlane(plane, 
										   data, 
										   axis,
										   &fGeom, 
										   &bGeom, 
										   pFront, 
										   pBack,
										   data.mIsKdNode);	 

	splitAxis = 3;

	if (candidateCost < lowestCost)
	{	
		bestPlane = plane;
		lowestCost = candidateCost;
		splitAxis = axis;
		// assign already computed values
		// we can do this because we always save the
		// computed values from the axis aligned splits		
		if (fGeom && bGeom)
		{
			frontData.mGeometry = fGeom;
			backData.mGeometry = bGeom;
	
			frontData.mProbability = pFront;
			backData.mProbability = pBack;
		}
	}
	else
	{
		DEL_PTR(fGeom);
		DEL_PTR(bGeom);
	}

	frontData.mIsKdNode = backData.mIsKdNode = 
		(data.mIsKdNode && splitAxis < 3);

#ifdef _DEBUG
	Debug << "plane lowest cost: " << lowestCost << endl;
#endif

	// cost ratio miss
	if (lowestCost > mTermMaxCostRatio)
	{
		return false;
	}

	return true;
}


Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
{
	const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));

	const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
	const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();

	const Vector3 pt = (maxPt + minPt) * 0.5;
	const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());

	return Plane3(normal, pt);
}


Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
{
	Vector3 pt[3];

	int idx[3];
	int cmaxT = 0;
	int cminT = 0;
	bool chooseMin = false;

	for (int j = 0; j < 3; ++ j)
	{
		idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));

		if (idx[j] >= (int)rays.size())
		{
			idx[j] -= (int)rays.size();

			chooseMin = (cminT < 2);
		}
		else
			chooseMin = (cmaxT < 2);

		RayInfo rayInf = rays[idx[j]];
		pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
	}

	return Plane3(pt[0], pt[1], pt[2]);
}


Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
{
	Vector3 pt[3];

	int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
	int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));

	// check if rays different
	if (idx1 == idx2)
		idx2 = (idx2 + 1) % (int)rays.size();

	const RayInfo ray1 = rays[idx1];
	const RayInfo ray2 = rays[idx2];

	// normal vector of the plane parallel to both lines
	const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));

	// vector from line 1 to line 2
	const Vector3 vd = ray2.ExtrapOrigin() - ray1.ExtrapOrigin();

	// project vector on normal to get distance
	const float dist = DotProd(vd, norm);

	// point on plane lies halfway between the two planes
	const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;

	return Plane3(norm, planePt);
}


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


float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
									 const VspBspTraversalData &data,
									 BspNodeGeometry &geomFront,
									 BspNodeGeometry &geomBack,
									 float &pFront,
									 float &pBack) const
{
	float cost = 0;

	float sumBalancedRays = 0;
	float sumRaySplits = 0;

	int pvsFront = 0;
	int pvsBack = 0;

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

	int raysFront = 0;
	int raysBack = 0;
	int totalPvs = 0;

	int limit;
	bool useRand;

	// choose test rays randomly if too much
	if ((int)data.mRays->size() > mMaxTests)
	{
		useRand = true;
		limit = mMaxTests;
	}
	else
	{
		useRand = false;
		limit = (int)data.mRays->size();
	}
	
	for (int i = 0; i < limit; ++ i)
	{
		const int testIdx = useRand ? 
			(int)RandomValue(0, (Real)((int)data.mRays->size() - 1)) : i;
		RayInfo rayInf = (*data.mRays)[testIdx];

		float t;
		VssRay *ray = rayInf.mRay;
		const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);

        if (0)
		{
			if (cf >= 0)
				++ raysFront;
			if (cf <= 0)
				++ raysBack;
		}

		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
		{
			sumBalancedRays += cf;
		}

		if (mSplitPlaneStrategy & BALANCED_RAYS)
		{
			if (cf == 0)
				++ sumRaySplits;
		}

		if (mSplitPlaneStrategy & PVS)
		{
			// find front and back pvs for origing and termination object
			AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
			AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
		}
	}

	const float raysSize = (float)data.mRays->size() + Limits::Small;

	if (mSplitPlaneStrategy & PVS)
	{
		// create unique ids for pvs heuristics
		GenerateUniqueIdsForPvs();

		// construct child geometry with regard to the candidate split plane
		data.mGeometry->SplitGeometry(geomFront,
									  geomBack,
									  candidatePlane,
									  mBox,
									  mEpsilon);

		pOverall = data.mProbability;

		if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
		{
			pFront = geomFront.GetVolume();
			pBack = pOverall - pFront;
		}
		else
		{
			pFront = geomFront.GetArea();
			pBack = geomBack.GetArea();
		}
	}


	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
		cost += mLeastRaySplitsFactor * sumRaySplits / raysSize;

	if (mSplitPlaneStrategy & BALANCED_RAYS)
		cost += mBalancedRaysFactor * fabs(sumBalancedRays) / raysSize;

	// -- pvs rendering heuristics
	if (mSplitPlaneStrategy & PVS)
	{
		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;

		float oldCost, newCost;

		// only render cost
		if (1)
		{
			oldCost = oldRenderCost;
			newCost = newRenderCost;
		}
		else // also considering standard deviation
		{
			// standard deviation is difference of back and front pvs
			const float expectedCost = 0.5f * (penaltyFront + penaltyBack);

			const float newDeviation = 0.5f * 
				fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);

			const float oldDeviation = penaltyOld;

			newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
			oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
		}

		cost += mPvsFactor * newCost / (oldCost + Limits::Small);
		
	}

#ifdef _DEBUG
  	Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
  		  << " frontpvs: " << pvsFront << " pFront: " << pFront
  		  << " backpvs: " << pvsBack << " pBack: " << pBack << endl << endl;
	Debug << "cost: " << cost << endl;
#endif

	
	// normalize cost by sum of linear factors
	if(0)
		return cost / (float)mCostNormalizer;
	else
		return cost;
}


float VspBspTree::EvalAxisAlignedSplitCost(const VspBspTraversalData &data,
										   const AxisAlignedBox3 &box,
										   const int axis,
										   const float &position,										  
										   float &pFront,
										   float &pBack) const
{
	int pvsTotal = 0;
	int pvsFront = 0;
	int 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)
	{
		//if (!(*rit).mRay->IsActive())	continue;

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

	if (!mUseAreaForPvs)
	{   // volume
		pBack = pFront = pOverall * 0.5f;
		
#if 0
		// box length substitute for probability
		const float minBox = box.Min(axis);
		const float maxBox = box.Max(axis);

		pBack = position - minBox;
		pFront = maxBox - position;
		pOverall = maxBox - minBox;
#endif
	}
	else //-- area substitute for probability
	{
		const int axis2 = (axis + 1) % 3;
		const int axis3 = (axis + 2) % 3;

		const float faceArea = 
			(box.Max(axis2) - box.Min(axis2)) *
			(box.Max(axis3) - box.Min(axis3));

		pBack = pFront = pOverall * 0.5f + faceArea;
	}

#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 VspBspTree::AddObjToPvs(Intersectable *obj,
							 const int cf,
							 int &frontPvs,
							 int &backPvs,
							 int &totalPvs) const
{
	if (!obj)
		return;
	
	if ((obj->mMailbox != sFrontId) &&
		(obj->mMailbox != sBackId) &&
		(obj->mMailbox != sFrontAndBackId))
	{
		++ totalPvs;
	}

	// 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;
		
			if (obj->mMailbox == sBackId)
				obj->mMailbox = sFrontAndBackId;
			else
				obj->mMailbox = sFrontId;
		}
	}

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


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

	while (!nodeStack.empty())
	{
		BspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
		{
			// test if this leaf is in valid view space
			BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
			if (leaf->TreeValid() && 
				(!onlyUnmailed || !leaf->Mailed()) &&
				((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
			{
				leaves.push_back(leaf);
			}
		}
		else
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

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


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


BspNode *VspBspTree::GetRoot() const
{
	return mRoot;
}


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

	// store maximal and minimal depth
	if (data.mDepth > mBspStats.maxDepth)
		mBspStats.maxDepth = data.mDepth;

	if (data.mPvs > mBspStats.maxPvs)
		mBspStats.maxPvs = data.mPvs;
	
	if (data.mDepth < mBspStats.minDepth)
		mBspStats.minDepth = data.mDepth;

	if (data.mDepth >= mTermMaxDepth)
		++ mBspStats.maxDepthNodes;

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

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

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

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

	if (data.mProbability <= mTermMinProbability)
		++ mBspStats.minProbabilityNodes;
	
	// accumulate depth to compute average depth
	mBspStats.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
}


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

	stack<BspRayTraversalData> tQueue;

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

	BspNode *node = mRoot;
	BspNode *farChild = NULL;

	while (1)
	{
		if (!node->IsLeaf())
		{
			BspInterior *in = dynamic_cast<BspInterior *>(node);

			Plane3 splitPlane = in->GetPlane();
			const int entSide = splitPlane.Side(entp);
			const int extSide = splitPlane.Side(extp);

			if (entSide < 0)
			{
				node = in->GetBack();

				if(extSide <= 0) // plane does not split ray => no far child
					continue;

				farChild = in->GetFront(); // plane splits ray

			} else if (entSide > 0)
			{
				node = in->GetFront();

				if (extSide >= 0) // plane does not split ray => no far child
					continue;

				farChild = in->GetBack(); // plane splits ray
			}
			else // ray and plane are coincident
			{
				// WHAT TO DO IN THIS CASE ?
				//break;
				node = in->GetFront();
				continue;
			}

			// push data for far child
			tQueue.push(BspRayTraversalData(farChild, extp, maxt));

			// find intersection of ray segment with plane
			float t;
			extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
			maxt *= t;

		} else // reached leaf => intersection with view cell
		{
			BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);

			if (!leaf->GetViewCell()->Mailed())
			{
				//ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
				leaf->GetViewCell()->Mail();
				++ hits;
			}

   			//-- fetch the next far child from the stack
			if (tQueue.empty())
				break;

			entp = extp;
			mint = maxt; // NOTE: need this?

			if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
				break;

			BspRayTraversalData &s = tQueue.top();

			node = s.mNode;
			extp = s.mExitPoint;
			maxt = s.mMaxT;

			tQueue.pop();
		}
	}

	return hits;
}


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


void VspBspTree::CollapseViewCells()
{
// TODO
#if VC_HISTORY
	stack<BspNode *> nodeStack;

	if (!mRoot)
		return;

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

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

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

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

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

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


void VspBspTree::CollectRays(VssRayContainer &rays)
{
	vector<BspLeaf *> leaves;

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

	for (lit = leaves.begin(); lit != lit_end; ++ lit)
	{
		BspLeaf *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 VspBspTree::ValidateTree()
{
	stack<BspNode *> nodeStack;

	if (!mRoot)
		return;

	nodeStack.push(mRoot);
	
	mBspStats.invalidLeaves = 0;
	while (!nodeStack.empty()) 
	{
		BspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
		{
			BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);

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

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

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



void VspBspTree::CollectViewCells(BspNode *root,
								  bool onlyValid,
								  ViewCellContainer &viewCells,
								  bool onlyUnmailed) const
{
	stack<BspNode *> nodeStack;

	if (!root)
		return;

	nodeStack.push(root);
	
	while (!nodeStack.empty()) 
	{
		BspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf())
		{
			if (!onlyValid || node->TreeValid())
			{
				ViewCell *viewCell = 
					mViewCellsTree->GetActiveViewCell(dynamic_cast<BspLeaf *>(node)->GetViewCell());
						
				if (!onlyUnmailed || !viewCell->Mailed()) 
				{
					viewCell->Mail();
					viewCells.push_back(viewCell);
				}
			}
		}
		else
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);
		
			nodeStack.push(interior->GetFront());
			nodeStack.push(interior->GetBack());
		}
	}

}


float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
{
	float len = 0;

	RayInfoContainer::const_iterator it, it_end = rays.end();

	for (it = rays.begin(); it != it_end; ++ it)
		len += (*it).SegmentLength();

	return len;
}


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

	RayInfoContainer::const_iterator it, it_end = rays.end();

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

		// get classification and receive new t
		const int cf = bRay.ComputeRayIntersection(plane, t);

		switch (cf)
		{
		case -1:
			backRays.push_back(bRay);
			break;
		case 1:
			frontRays.push_back(bRay);
			break;
		case 0:
			{
				//-- split ray
				//   test if start point behind or in front of plane
				const int side = plane.Side(bRay.ExtrapOrigin());

				if (side <= 0)
				{
					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()));
				}
			}
			break;
		default:
			Debug << "Should not come here" << endl;
			break;
		}
	}

	return splits;
}


void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
{
	BspNode *lastNode;

	do
	{
		lastNode = n;

		// want to get planes defining geometry of this node => don't take
		// split plane of node itself
		n = n->GetParent();

		if (n)
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(n);
			Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();

            if (interior->GetFront() != lastNode)
				halfSpace.ReverseOrientation();

			halfSpaces.push_back(halfSpace);
		}
	}
	while (n);
}


void VspBspTree::ConstructGeometry(BspNode *n,
								   BspNodeGeometry &geom) const
{
	vector<Plane3> halfSpaces;
	ExtractHalfSpaces(n, halfSpaces);

	PolygonContainer candidatePolys;

	// bounded planes are added to the polygons (reverse polygons
	// as they have to be outfacing
	for (int i = 0; i < (int)halfSpaces.size(); ++ i)
	{
		Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);

		if (p->Valid(mEpsilon))
		{
			candidatePolys.push_back(p->CreateReversePolygon());
			DEL_PTR(p);
		}
	}

	// add faces of bounding box (also could be faces of the cell)
	for (int i = 0; i < 6; ++ i)
	{
		VertexContainer vertices;

		for (int j = 0; j < 4; ++ j)
			vertices.push_back(mBox.GetFace(i).mVertices[j]);

		candidatePolys.push_back(new Polygon3(vertices));
	}

	for (int i = 0; i < (int)candidatePolys.size(); ++ i)
	{
		// polygon is split by all other planes
		for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
		{
			if (i == j) // polygon and plane are coincident
				continue;

			VertexContainer splitPts;
			Polygon3 *frontPoly, *backPoly;

			const int cf =
				candidatePolys[i]->ClassifyPlane(halfSpaces[j],
												 mEpsilon);

			switch (cf)
			{
				case Polygon3::SPLIT:
					frontPoly = new Polygon3();
					backPoly = new Polygon3();

					candidatePolys[i]->Split(halfSpaces[j],
											 *frontPoly,
											 *backPoly,
											 mEpsilon);

					DEL_PTR(candidatePolys[i]);

					if (frontPoly->Valid(mEpsilon))
						candidatePolys[i] = frontPoly;
					else
						DEL_PTR(frontPoly);

					DEL_PTR(backPoly);
					break;
				case Polygon3::BACK_SIDE:
					DEL_PTR(candidatePolys[i]);
					break;
				// just take polygon as it is
				case Polygon3::FRONT_SIDE:
				case Polygon3::COINCIDENT:
				default:
					break;
			}
		}

		if (candidatePolys[i])
			geom.mPolys.push_back(candidatePolys[i]);
	}
}


void VspBspTree::ConstructGeometry(ViewCell *vc, 
								   BspNodeGeometry &vcGeom) const
{
	ViewCellContainer leaves;
	
	mViewCellsTree->CollectLeaves(vc, leaves);

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

	for (it = leaves.begin(); it != it_end; ++ it)
	{
		BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
		
		ConstructGeometry(l, vcGeom);
	}
}


typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;


int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
							  const bool onlyUnmailed) const
{
	stack<bspNodePair> nodeStack;
	
	BspNodeGeometry nodeGeom;
	ConstructGeometry(n, nodeGeom);
	
	// split planes from the root to this node
	// needed to verify that we found neighbor leaf
	// TODO: really needed?
	vector<Plane3> halfSpaces;
	ExtractHalfSpaces(n, halfSpaces);


	BspNodeGeometry *rgeom = new BspNodeGeometry();
	ConstructGeometry(mRoot, *rgeom);

	nodeStack.push(bspNodePair(mRoot, rgeom));

	while (!nodeStack.empty())
	{
		BspNode *node = nodeStack.top().first;
		BspNodeGeometry *geom = nodeStack.top().second;
	
		nodeStack.pop();

		if (node->IsLeaf())
		{
			// test if this leaf is in valid view space
			if (node->TreeValid() &&
				(node != n) && 
				(!onlyUnmailed || !node->Mailed()))
			{
				bool isAdjacent = true;

				if (1)
				{
					// test all planes of current node if still adjacent
					for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
					{
						const int cf =
							Polygon3::ClassifyPlane(geom->mPolys,
													halfSpaces[i],
													mEpsilon);

						if (cf == Polygon3::BACK_SIDE)
						{
							isAdjacent = false;
						}
					}
				}
				else
				{
					// TODO: why is this wrong??
					// test all planes of current node if still adjacent
					for (int i = 0; (i < (int)nodeGeom.mPolys.size()) && isAdjacent; ++ i)
					{
						Polygon3 *poly = nodeGeom.mPolys[i];

						const int cf =
							Polygon3::ClassifyPlane(geom->mPolys,
													poly->GetSupportingPlane(),
													mEpsilon);

						if (cf == Polygon3::BACK_SIDE)
						{
							isAdjacent = false;
						}
					}
				}
				// neighbor was found
				if (isAdjacent)
				{	
					neighbors.push_back(dynamic_cast<BspLeaf *>(node));
				}
			}
		}
		else
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

			const int cf = Polygon3::ClassifyPlane(nodeGeom.mPolys,
												   interior->GetPlane(),
												   mEpsilon);
			
			BspNode *front = interior->GetFront();
			BspNode *back = interior->GetBack();
            
			BspNodeGeometry *fGeom = new BspNodeGeometry();
			BspNodeGeometry *bGeom = new BspNodeGeometry();

			geom->SplitGeometry(*fGeom,
								*bGeom,
								interior->GetPlane(),
								mBox,
								mEpsilon);
		
			if (cf == Polygon3::FRONT_SIDE)
			{
				nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
				DEL_PTR(bGeom);
			}
			else
			{
				if (cf == Polygon3::BACK_SIDE)
				{
					nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
					DEL_PTR(fGeom);
				}
				else
				{	// random decision
					nodeStack.push(bspNodePair(front, fGeom));
					nodeStack.push(bspNodePair(back, bGeom));
				}
			}
		}
	
		DEL_PTR(geom);
	}

	return (int)neighbors.size();
}



int VspBspTree::FindApproximateNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
										 const bool onlyUnmailed) const
{
	stack<bspNodePair> nodeStack;
	
	BspNodeGeometry nodeGeom;
	ConstructGeometry(n, nodeGeom);
	
	// split planes from the root to this node
	// needed to verify that we found neighbor leaf
	// TODO: really needed?
	vector<Plane3> halfSpaces;
	ExtractHalfSpaces(n, halfSpaces);


	BspNodeGeometry *rgeom = new BspNodeGeometry();
	ConstructGeometry(mRoot, *rgeom);

	nodeStack.push(bspNodePair(mRoot, rgeom));

	while (!nodeStack.empty())
	{
		BspNode *node = nodeStack.top().first;
		BspNodeGeometry *geom = nodeStack.top().second;
	
		nodeStack.pop();

		if (node->IsLeaf())
		{
			// test if this leaf is in valid view space
			if (node->TreeValid() &&
				(node != n) && 
				(!onlyUnmailed || !node->Mailed()))
			{
				bool isAdjacent = true;

				// neighbor was found
				if (isAdjacent)
				{	
					neighbors.push_back(dynamic_cast<BspLeaf *>(node));
				}
			}
		}
		else
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

			const int cf = Polygon3::ClassifyPlane(nodeGeom.mPolys,
												   interior->GetPlane(),
												   mEpsilon);
			
			BspNode *front = interior->GetFront();
			BspNode *back = interior->GetBack();
            
			BspNodeGeometry *fGeom = new BspNodeGeometry();
			BspNodeGeometry *bGeom = new BspNodeGeometry();

			geom->SplitGeometry(*fGeom,
								*bGeom,
								interior->GetPlane(),
								mBox,
								mEpsilon);
		
			if (cf == Polygon3::FRONT_SIDE)
			{
				nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
				DEL_PTR(bGeom);
			}
			else
			{
				if (cf == Polygon3::BACK_SIDE)
				{
					nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
					DEL_PTR(fGeom);
				}
				else
				{	// random decision
					nodeStack.push(bspNodePair(front, fGeom));
					nodeStack.push(bspNodePair(back, bGeom));
				}
			}
		}
	
		DEL_PTR(geom);
	}

	return (int)neighbors.size();
}



BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
{
    stack<BspNode *> nodeStack;
  	nodeStack.push(mRoot);

	int mask = rand();

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

		if (node->IsLeaf())
		{
			return dynamic_cast<BspLeaf *>(node);
		}
		else
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);
			BspNode *next;
			BspNodeGeometry geom;

			// todo: not very efficient: constructs full cell everytime
			ConstructGeometry(interior, geom);

			const int cf = 
				Polygon3::ClassifyPlane(geom.mPolys, halfspace, mEpsilon);

			if (cf == Polygon3::BACK_SIDE)
				next = interior->GetFront();
			else
				if (cf == Polygon3::FRONT_SIDE)
					next = interior->GetFront();
			else
			{
				// random decision
				if (mask & 1)
					next = interior->GetBack();
				else
					next = interior->GetFront();
				mask = mask >> 1;
			}

			nodeStack.push(next);
		}
	}

	return NULL;
}

BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
{
	stack<BspNode *> nodeStack;

	nodeStack.push(mRoot);

	int mask = rand();

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

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

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

			mask = mask >> 1;
		}
	}

	return NULL;
}

int VspBspTree::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 VspBspTree::GetEpsilon() const
{
	return mEpsilon;
}


int VspBspTree::SplitPolygons(const Plane3 &plane,
							  PolygonContainer &polys,
							  PolygonContainer &frontPolys,
							  PolygonContainer &backPolys,
							  PolygonContainer &coincident) const
{
	int splits = 0;

	PolygonContainer::const_iterator it, it_end = polys.end();

	for (it = polys.begin(); it != polys.end(); ++ it)	
	{
		Polygon3 *poly = *it;

		// classify polygon
		const int cf = poly->ClassifyPlane(plane, mEpsilon);

		switch (cf)
		{
			case Polygon3::COINCIDENT:
				coincident.push_back(poly);
				break;
			case Polygon3::FRONT_SIDE:
				frontPolys.push_back(poly);
				break;
			case Polygon3::BACK_SIDE:
				backPolys.push_back(poly);
				break;
			case Polygon3::SPLIT:
				backPolys.push_back(poly);
				frontPolys.push_back(poly);
				++ splits;
				break;
			default:
                Debug << "SHOULD NEVER COME HERE\n";
				break;
		}
	}

	return splits;
}


int VspBspTree::CastLineSegment(const Vector3 &origin,
								const Vector3 &termination,
								vector<ViewCell *> &viewcells)
{
	int hits = 0;
	stack<BspRayTraversalData> tQueue;

	float mint = 0.0f, maxt = 1.0f;

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

	Vector3 entp = origin;
	Vector3 extp = termination;

	BspNode *node = mRoot;
	BspNode *farChild = NULL;

	float t;
	while (1)
	{
		if (!node->IsLeaf())
		{
			BspInterior *in = dynamic_cast<BspInterior *>(node);

			Plane3 splitPlane = in->GetPlane();
			
			const int entSide = splitPlane.Side(entp);
			const int extSide = splitPlane.Side(extp);

			if (entSide < 0) 
			{
			  node = in->GetBack();
			  // plane does not split ray => no far child
			  if (extSide <= 0) 
				continue;
			  
			  farChild = in->GetFront(); // plane splits ray
			} 
			else if (entSide > 0)
			{
			  node = in->GetFront();

			  if (extSide >= 0) // plane does not split ray => no far child
				continue;

				farChild = in->GetBack(); // plane splits ray
			}
			else // ray end point on plane
			{	// NOTE: what to do if ray is coincident with plane?
				if (extSide < 0)
					node = in->GetBack();
				else 
					node = in->GetFront();
								
				continue; // no far child
			}

			// push data for far child
			tQueue.push(BspRayTraversalData(farChild, extp));

			// find intersection of ray segment with plane
			extp = splitPlane.FindIntersection(origin, extp, &t);
		} 
		else
		{
			// reached leaf => intersection with view cell
			BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);

			ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
			if (!viewCell->Mailed())
			{
				viewcells.push_back(viewCell);
				viewCell->Mail();
				++ hits;
			}

			//-- fetch the next far child from the stack
			if (tQueue.empty())
				break;

			entp = extp;
			
			BspRayTraversalData &s = tQueue.top();

			node = s.mNode;
			extp = s.mExitPoint;

			tQueue.pop();
		}
	}

	return hits;
}




int VspBspTree::TreeDistance(BspNode *n1, BspNode *n2) const
{
	std::deque<BspNode *> path1;
	BspNode *p1 = n1;

	// create path from node 1 to root
	while (p1)
	{
		if (p1 == n2) // second node on path
			return (int)path1.size();

		path1.push_front(p1);
		p1 = p1->GetParent();
	}

	int depth = n2->GetDepth();
	int d = depth;

	BspNode *p2 = n2;

	// compare with same depth
	while (1)
	{
		if ((d < (int)path1.size()) && (p2 == path1[d]))
			return (depth - d) + ((int)path1.size() - 1 - d);

		-- d;
		p2 = p2->GetParent();
	}

	return 0; // never come here
}


BspNode *VspBspTree::CollapseTree(BspNode *node, int &collapsed)
{
// TODO
#if VC_HISTORY
	if (node->IsLeaf())
		return node;

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

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

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

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

			BspLeaf *leaf = new BspLeaf(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 VspBspTree::CollapseTree()
{
	int collapsed = 0;
	//TODO
#if VC_HISTORY
	(void)CollapseTree(mRoot, collapsed);

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


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

	ViewCell::NewMail();

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

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

			BspViewCell *viewCell = leaf->GetViewCell();

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

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

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


typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;


int VspBspTree::CastBeam(Beam &beam)
{
    stack<bspNodePair> nodeStack;
	BspNodeGeometry *rgeom = new BspNodeGeometry();
	ConstructGeometry(mRoot, *rgeom);

	nodeStack.push(bspNodePair(mRoot, rgeom));
  
	ViewCell::NewMail();

	while (!nodeStack.empty()) 
	{
		BspNode *node = nodeStack.top().first;
		BspNodeGeometry *geom = nodeStack.top().second;
		nodeStack.pop();
		
		AxisAlignedBox3 box;
		box.Initialize();
		geom->IncludeInBox(box);

		const int side = beam.ComputeIntersection(box);
		
		switch (side) 
		{
		case -1:
			CollectViewCells(node, true, beam.mViewCells, true);
			break;
		case 0:
			
			if (node->IsLeaf())
			{
				BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
			
				if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
				{
					leaf->GetViewCell()->Mail();
					beam.mViewCells.push_back(leaf->GetViewCell());
				}
			}
			else 
			{
				BspInterior *interior = dynamic_cast<BspInterior *>(node);
			
				BspNode *first = interior->GetFront();
				BspNode *second = interior->GetBack();
            
				BspNodeGeometry *firstGeom = new BspNodeGeometry();
				BspNodeGeometry *secondGeom = new BspNodeGeometry();

				geom->SplitGeometry(*firstGeom,
									*secondGeom,
									interior->GetPlane(),
									mBox,
									mEpsilon);

				// decide on the order of the nodes
				if (DotProd(beam.mPlanes[0].mNormal, 
					interior->GetPlane().mNormal) > 0)
				{
					swap(first, second);
					swap(firstGeom, secondGeom);
				}

				nodeStack.push(bspNodePair(first, firstGeom));
				nodeStack.push(bspNodePair(second, secondGeom));
			}
			
			break;
		default:
			// default: cull
			break;
		}
		
		DEL_PTR(geom);
		
	}

	return (int)beam.mViewCells.size();
}


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


int VspBspTree::CollectMergeCandidates(const vector<BspLeaf *> leaves,
									   vector<MergeCandidate> &candidates)
{
	BspLeaf::NewMail();
	
	vector<BspLeaf *>::const_iterator it, it_end = leaves.end();

	int numCandidates = 0;

	// find merge candidates and push them into queue
	for (it = leaves.begin(); it != it_end; ++ it)
	{
		BspLeaf *leaf = *it;
		
		// the same leaves must not be part of two merge candidates
		leaf->Mail();
		vector<BspLeaf *> neighbors;
		FindNeighbors(leaf, neighbors, true);
		//FindApproximateNeighbors(leaf, neighbors, true);
		vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();

		// TODO: test if at least one ray goes from one leaf to the other
		for (nit = neighbors.begin(); nit != nit_end; ++ nit)
		{
			if ((*nit)->GetViewCell() != leaf->GetViewCell())
			{
				MergeCandidate mc(leaf->GetViewCell(), (*nit)->GetViewCell());
				candidates.push_back(mc);

				++ numCandidates;
				if ((numCandidates % 1000) == 0)
				{
					cout << "collected " << numCandidates << " merge candidates" << endl;
				}
			}
		}
	}

	Debug << "merge queue: " << (int)candidates.size() << endl;
	Debug << "leaves in queue: " << numCandidates << endl;
	

	return (int)leaves.size();
}


int VspBspTree::CollectMergeCandidates(const VssRayContainer &rays, 
									   vector<MergeCandidate> &candidates)
{
	ViewCell::NewMail();
	long startTime = GetTime();
	
	map<BspLeaf *, vector<BspLeaf*> > neighborMap;
	ViewCellContainer::const_iterator iit;

	int numLeaves = 0;
	
	BspLeaf::NewMail();

	for (int i = 0; i < (int)rays.size(); ++ i)
	{  
		VssRay *ray = rays[i];
	
		// traverse leaves stored in the rays and compare and 
		// merge consecutive leaves (i.e., the neighbors in the tree)
		if (ray->mViewCells.size() < 2)
			continue;
//TODO viewcellhierarchy
		iit = ray->mViewCells.begin();
		BspViewCell *bspVc = dynamic_cast<BspViewCell *>(*(iit ++));
		BspLeaf *leaf = bspVc->mLeaf;
		
		// traverse intersections 
		// consecutive leaves are neighbors => add them to queue
		for (; iit != ray->mViewCells.end(); ++ iit)
		{
			// next pair
			BspLeaf *prevLeaf = leaf;
			bspVc = dynamic_cast<BspViewCell *>(*iit);
            leaf = bspVc->mLeaf;

			// view space not valid or same view cell
			if (!leaf->TreeValid() || !prevLeaf->TreeValid() ||
				(leaf->GetViewCell() == prevLeaf->GetViewCell()))
				continue;

           	vector<BspLeaf *> &neighbors = neighborMap[leaf];
			
			bool found = false;

			// both leaves inserted in queue already =>
			// look if double pair already exists
			if (leaf->Mailed() && prevLeaf->Mailed())
			{
				vector<BspLeaf *>::const_iterator it, it_end = neighbors.end();
				
                for (it = neighbors.begin(); !found && (it != it_end); ++ it)
					if (*it == prevLeaf)
						found = true; // already in queue
			}
		
			if (!found)
			{
				// this pair is not in map yet
				// => insert into the neighbor map and the queue
				neighbors.push_back(prevLeaf);
				neighborMap[prevLeaf].push_back(leaf);

				leaf->Mail();
				prevLeaf->Mail();
		
				MergeCandidate mc(leaf->GetViewCell(), prevLeaf->GetViewCell());
				
				candidates.push_back(mc);

				if (((int)candidates.size() % 1000) == 0)
				{
					cout << "collected " << (int)candidates.size() << " merge candidates" << endl;
				}
			}
        }
	}

	Debug << "neighbormap size: " << (int)neighborMap.size() << endl;
	Debug << "merge queue: " << (int)candidates.size() << endl;
	Debug << "leaves in queue: " << numLeaves << endl;


	//-- collect the leaves which haven't been found by ray casting
	if (0)
	{
		cout << "finding additional merge candidates using geometry" << endl;
		vector<BspLeaf *> leaves;
		CollectLeaves(leaves, true);
		Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
		CollectMergeCandidates(leaves, candidates);
	}

	return numLeaves;
}




ViewCell *VspBspTree::GetViewCell(const Vector3 &point)
{
  if (mRoot == NULL)
	return NULL;
  
  stack<BspNode *> nodeStack;
  nodeStack.push(mRoot);
  
  ViewCell *viewcell = NULL;
  
  while (!nodeStack.empty())  {
	BspNode *node = nodeStack.top();
	nodeStack.pop();
	
	if (node->IsLeaf()) {
	  viewcell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
	  break;
	} else {
	  
	  BspInterior *interior = dynamic_cast<BspInterior *>(node);
		
	  // random decision
	  if (interior->GetPlane().Side(point) < 0)
		nodeStack.push(interior->GetBack());
	  else
		nodeStack.push(interior->GetFront());
	}
  }
  
  return viewcell;
}


bool VspBspTree::ViewPointValid(const Vector3 &viewPoint) const
{
	BspNode *node = mRoot;

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

		if (node->IsLeaf())
			return false;
			
		BspInterior *in = dynamic_cast<BspInterior *>(node);
					
		if (in->GetPlane().Side(viewPoint) <= 0) 
		{
			node = in->GetBack();
		} 
		else
		{
			node = in->GetFront();
		}
	}

	// should never come here
	return false;
}


void VspBspTree::PropagateUpValidity(BspNode *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();
			BspInterior *interior = dynamic_cast<BspInterior *>(node);
			
			// the parent is valid iff both leaves are valid
			node->SetTreeValid(interior->GetBack()->TreeValid() && 
							   interior->GetFront()->TreeValid());
		}
	}
}


bool VspBspTree::Export(ofstream &stream)
{
	ExportNode(mRoot, stream);

	return true;
}


void VspBspTree::ExportNode(BspNode *node, ofstream &stream)
{
	if (node->IsLeaf())
	{
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
		ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());

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

		stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
	}
	else
	{
		BspInterior *interior = dynamic_cast<BspInterior *>(node);
	
		Plane3 plane = interior->GetPlane();
		stream << "<Interior plane=\"" << plane.mNormal.x << " " 
			   << plane.mNormal.y << " " << plane.mNormal.z << " " 
			   << plane.mD << "\">" << endl;

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

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