#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"
#include "IntersectableWrapper.h"



namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0
#define COUNT_ORIGIN_OBJECTS 1


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

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



typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;


// pvs penalty can be different from pvs size
inline static float EvalPvsPenalty(const float pvs, 
								   const float lower,
								   const float upper)
{
	// clamp to minmax values
	if (pvs < lower)
		return (float)lower;
	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(1)
{
	bool randomize = false;
	Environment::GetSingleton()->GetBoolValue("VspBspTree.Construction.randomize", randomize);
	if (randomize) Randomize(); // initialise random generator for heuristics

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

	Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
	Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.minProbability", mTermMinProbability);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
	
	Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
	Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);

	////////////////////////
	//-- cost ratios for early tree termination
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
	Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);

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

	Environment::GetSingleton()->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);

	//////////
	//-- partition criteria

	Environment::GetSingleton()->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
	Environment::GetSingleton()->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
	Environment::GetSingleton()->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);

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

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

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

	Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.renderCostWeight", mRenderCostWeight);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
	Environment::GetSingleton()->GetBoolValue("VspBspTree.usePolygonSplitIfAvailable", mUsePolygonSplitIfAvailable);

	Environment::GetSingleton()->GetBoolValue("VspBspTree.useCostHeuristics", mUseCostHeuristics);
	Environment::GetSingleton()->GetBoolValue("VspBspTree.useSplitCostQueue", mUseSplitCostQueue);
	Environment::GetSingleton()->GetBoolValue("VspBspTree.simulateOctree", mCirculatingAxis);
	Environment::GetSingleton()->GetBoolValue("VspBspTree.useRandomAxis", mUseRandomAxis);
	Environment::GetSingleton()->GetIntValue("VspBspTree.nodePriorityQueueType", mNodePriorityQueueType);

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

	Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.minBand", mMinBand);
	Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.maxBand", mMaxBand);
	Environment::GetSingleton()->GetBoolValue("VspBspTree.Construction.useDrivingAxisForMaxCost", mUseDrivingAxisIfMaxCostViolated);

	/////////
	//-- 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 << "max polygon candidates: " << mMaxPolyCandidates << endl;
	Debug << "randomize: " << randomize << endl;

	Debug << "using area for pvs: " << mUseAreaForPvs << endl;
	Debug << "render cost weight: " << mRenderCostWeight << 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 poly split if available: " << mUsePolygonSplitIfAvailable << endl;
	Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
	Debug << "use split cost queue: " << mUseSplitCostQueue << 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;
	Debug << "use driving axis for max cost: " << mUseDrivingAxisIfMaxCostViolated << endl;
	Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << 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";
	}

	Debug << endl;

	mLocalSubdivisionCandidates = new vector<SortableEntry>;
}


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


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

	return mOutOfBoundsCell;
}


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


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


int VspBspTree::AddMeshToPolygons(Mesh *mesh,
								  PolygonContainer &polys) const
{
	if (!mesh) return 0;

	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))
		{
			polys.push_back(poly);
		}
		else
		{
			DEL_PTR(poly);
		}
	}

	return (int)mesh->mFaces.size();
}


void VspBspTree::ExtractPolygons(Intersectable *object, PolygonContainer &polys) const
{
	// extract the polygons from the intersectables
	switch (object->Type())
	{
		case Intersectable::MESH_INSTANCE:
			{
				Mesh *mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
				AddMeshToPolygons(mesh, polys);
			}
			break;
		case Intersectable::VIEW_CELL:
			{
				Mesh *mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
				AddMeshToPolygons(mesh, polys);
				break;
			}
		case Intersectable::TRANSFORMED_MESH_INSTANCE:
			{
				TransformedMeshInstance *mi = 
					dynamic_cast<TransformedMeshInstance *>(object);

				Mesh mesh;
				mi->GetTransformedMesh(mesh);
				AddMeshToPolygons(&mesh, polys);
			}
			break;
		case Intersectable::TRIANGLE_INTERSECTABLE:
			{
				TriangleIntersectable *intersect = 
					dynamic_cast<TriangleIntersectable *>(object);

				Polygon3 *poly = new Polygon3(intersect->GetItem());

				if (poly->Valid(mEpsilon))	
				{
					polys.push_back(poly);
				}
				else
				{
					delete poly;
				}
			}
			break;
		default:
			Debug << "intersectable type not supported" << endl;
			break;
	}
}


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

	for (int i = 0; i < limit; ++i)
	{
		Intersectable *object = objects[i];//*it;
		ExtractPolygons(object, polys);

		mBoundingBox.Include(object->GetBox()); // add to BSP tree aabb
	}

	return (int)polys.size();
}


void VspBspTree::ComputeBoundingBox(const VssRayContainer &sampleRays,
									AxisAlignedBox3 *forcedBoundingBox) 
{
	if (forcedBoundingBox)
	{
		mBoundingBox = *forcedBoundingBox;
	}
	else // compute vsp tree bounding box
	{
		mBoundingBox.Initialize();

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

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

			// compute bounding box of view space
			mBoundingBox.Include(ray->GetTermination());
			mBoundingBox.Include(ray->GetOrigin());
		}
	}
}


void VspBspTree::Construct(const VssRayContainer &sampleRays,
						   AxisAlignedBox3 *forcedBoundingBox)
{
	// Compute the bounding box from the rays
	ComputeBoundingBox(sampleRays, forcedBoundingBox);
	
	PolygonContainer polys;
	RayInfoContainer *rays = new RayInfoContainer();

	////////////
	//-- extract polygons from rays if polygon candidate planes are required

	if (mMaxPolyCandidates) 
	{
		int numObj = 0;
		Intersectable::NewMail();

        cout << "Extracting polygons from rays ... ";
		const long startTime = GetTime();

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

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

			++ numObj;

			/////////
			//-- compute bounding box

			if (!forcedBoundingBox)
			{
				mBoundingBox.Include(ray->mTermination);
			}

			if ((mBoundingBox.IsInside(ray->mOrigin) || !forcedBoundingBox) &&
				ray->mOriginObject && 
				!ray->mOriginObject->Mailed())
			{		
				ray->mOriginObject->Mail();
				ExtractPolygons(ray->mOriginObject, polys);
								
				++ numObj;
			}
		}

		// throw out unnecessary polygons
		PreprocessPolygons(polys);
		cout << "finished" << endl;

		Debug << "\n" << (int)polys.size() << " polys extracted from "
		  << (int)sampleRays.size() << " rays in "
		  << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;
	}
	
	Debug << "maximal pvs (i.e., pvs still considered as valid): " 
		  << mViewCellsManager->GetMaxPvsSize() << endl;

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

	//-- 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 (mBoundingBox.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 *= mBoundingBox.SurfaceArea(); 
	else 
		mTermMinProbability *= mBoundingBox.GetVolume();


	mBspStats.nodes = 1;
	mBspStats.polys = (int)polys.size();
	mBspStats.mGlobalCostMisses = 0;


	// use split cost priority queue
	if (mUseSplitCostQueue)
	{
		ConstructWithSplitQueue(polys, rays);
	}
	else
	{
		Construct(polys, rays);
	}

	// clean up polygons
	CLEAR_CONTAINER(polys);
}


// TODO: return memory usage in MB
float VspBspTree::GetMemUsage() const
{
	return (float)
		 (sizeof(VspBspTree) + 
		  mBspStats.Leaves() * sizeof(BspLeaf) + 
		  mCreatedViewCells * sizeof(BspViewCell) +
		  mBspStats.pvs * sizeof(PvsData) +
		  mBspStats.Interior() * sizeof(BspInterior) +
		  mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
}


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

	/// create new vsp tree
	mRoot = new BspLeaf();

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

	/// we use the overall probability as normalizer
	/// either the overall area or the volume
	const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();

	/// first traversal data
	VspBspTraversalData tData(mRoot,
							  new PolygonContainer(polys),
							  0,
							  rays,
                              ComputePvsSize(*rays),
							  prop,
							  geom);

	// evaluate the priority of this traversal data
	EvalPriority(tData);

	// 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 / mBoundingBox.GetVolume();
	mTotalPvsSize = tData.mPvs;
	
	// first subdivison statistics
	AddSubdivisionStats(1, 0, 0, mTotalCost, (float)mTotalPvsSize);
    
	mBspStats.Start();
	cout << "Constructing vsp bsp tree ... \n";

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

	int nLeaves = 500;
	int nViewCells = 500;

	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
		const 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 in " << TimeDiff(startTime, GetTime())*1e-3 << "secs" << endl;

	mBspStats.Stop();
}



void VspBspTree::ConstructWithSplitQueue(const PolygonContainer &polys, 
											  RayInfoContainer *rays)
{
	VspBspSplitQueue 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);


	// first node is kd node, i.e. an axis aligned box
	if (1)
        tData.mIsKdNode = true;
	else
		tData.mIsKdNode = false;

	// compute first split candidate
	VspBspSubdivisionCandidate splitCandidate;
	splitCandidate.mParentData = tData;

	EvalSubdivisionCandidate(splitCandidate);

	tQueue.push(splitCandidate);

	mTotalCost = tData.mPvs * tData.mProbability / mBoundingBox.GetVolume();
	mTotalPvsSize = tData.mPvs;
	
	// first subdivison statistics
	AddSubdivisionStats(1, 0, 0, mTotalCost, (float)mTotalPvsSize);
   
    mBspStats.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.mRenderCostDecr / mTotalCost;

		//Debug << "cost ratio: " << costRatio << endl;
		if (costRatio < mTermMinGlobalCostRatio)
		{
			++ mBspStats.mGlobalCostMisses;
		}

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

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

		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::LocalTerminationCriteriaMet(const VspBspTraversalData &data) const
{
	return
		(((int)data.mRays->size() <= mTermMinRays) ||
		 (data.mPvs <= mTermMinPvs)   ||
		 (data.mProbability <= mTermMinProbability) ||
		 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
		 (data.mDepth >= mTermMaxDepth));
}


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


bool VspBspTree::GlobalTerminationCriteriaMet(const VspBspTraversalData &data) const
{
	return
		(0
		|| mOutOfMemory 
		|| (mBspStats.Leaves() >= mMaxViewCells) 
		|| (mBspStats.mGlobalCostMisses >= mTermGlobalCostMissTolerance)
		 );
}


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

	if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
	{
		PolygonContainer coincident;

		VspBspTraversalData tFrontData;
		VspBspTraversalData tBackData;

		// create new interior node and two leaf nodes
		// or return leaf as it is (if maxCostRatio missed)
		int splitAxis;
		bool splitFurther = true;
		int maxCostMisses = tData.mMaxCostMisses;
		
		Plane3 splitPlane;
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
		
		// choose next split plane
		if (!SelectPlane(splitPlane, leaf, tData, tFrontData, tBackData, splitAxis))
		{
			++ maxCostMisses;

			if (maxCostMisses > mTermMissTolerance)
			{
				// terminate branch because of max cost
				++ mBspStats.maxCostNodes;
				splitFurther = false;
			}
		}
	
		// if this a valid split => subdivide this node further

		if (splitFurther)
		{
			newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);

			if (splitAxis < 3)
				++ mBspStats.splits[splitAxis];
			else
				++ mBspStats.polySplits;

			// if it was a kd node (i.e., a box) and the split axis is axis aligned, it is still a kd node
			tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && (splitAxis < 3));
			
			tFrontData.mAxis = tBackData.mAxis = splitAxis;

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

			EvalPriority(tFrontData);
			EvalPriority(tBackData);

			// evaluate subdivision stats
			if (1)
				EvalSubdivisionStats(tData, tFrontData, tBackData);
			

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

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

		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->mLeaves.push_back(leaf);

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

		leaf->mProbability = tData.mProbability;

		// finally evaluate stats until this leaf
		if (0)
			EvaluateLeafStats(tData);		
	}

	//-- cleanup
	tData.Clear();

	return newNode;
}


// subdivide node using a split plane queue
BspNode *VspBspTree::Subdivide(VspBspSplitQueue &tQueue,
							   VspBspSubdivisionCandidate &splitCandidate)
{
	VspBspTraversalData &tData = splitCandidate.mParentData;

	BspNode *newNode = tData.mNode;

	if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
	{	
		PolygonContainer coincident;

		VspBspTraversalData tFrontData;
		VspBspTraversalData tBackData;

		////////////////////
		//-- continue subdivision
		
		// create new interior node and two leaf node
		const Plane3 splitPlane = splitCandidate.mSplitPlane;
				
		newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
	
		const int splitAxis = splitCandidate.mSplitAxis;
		const int maxCostMisses = splitCandidate.mMaxCostMisses;

		if (splitAxis < 3)
			++ mBspStats.splits[splitAxis];
		else
			++ mBspStats.polySplits;

		tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && (splitAxis < 3));
		tFrontData.mAxis = tBackData.mAxis = splitAxis;

		// how often was max cost ratio missed in this branch?
		tFrontData.mMaxCostMisses = maxCostMisses;
		tBackData.mMaxCostMisses = maxCostMisses;
			
		// statistics
		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) / mBoundingBox.GetVolume();

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

			AddSubdivisionStats(mBspStats.Leaves(), 
								-costDecr,  
								splitCandidate.GetPriority(), 
								mTotalCost, 
								(float)mTotalPvsSize / (float)mBspStats.Leaves());
		}

	
		//-- push the new split candidates on the stack
		VspBspSubdivisionCandidate frontCandidate;
		frontCandidate.mParentData = tFrontData;

		VspBspSubdivisionCandidate backCandidate;
		backCandidate.mParentData = tBackData;

		EvalSubdivisionCandidate(frontCandidate);
		EvalSubdivisionCandidate(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())
	{
		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);

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

		mBspStats.contributingSamples += conSamp;
		mBspStats.sampleContributions += (int)sampCon;

		viewCell->mLeaves.push_back(leaf);

		///////////
		//-- 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);
			}
		}
	
		if (mUseAreaForPvs)
			viewCell->SetArea(tData.mProbability);
		else
			viewCell->SetVolume(tData.mProbability);

        leaf->mProbability = tData.mProbability;

		// finally evaluate stats until this leaf
		if (0)
			EvaluateLeafStats(tData);		
	}

	//-- cleanup
	tData.Clear();

	return newNode;
}


void VspBspTree::EvalPriority(VspBspTraversalData &tData) const
{
    switch (mNodePriorityQueueType) 
	{
	case BREATH_FIRST:
		tData.mPriority = (float)-tData.mDepth;
		break;
	case DEPTH_FIRST:
		tData.mPriority = (float)tData.mDepth;
		break;
	default:
		tData.mPriority = tData.mPvs * tData.mProbability;
		//Debug << "priority: " << tData.mPriority << endl;
		break;
	}
}


void VspBspTree::EvalSubdivisionCandidate(VspBspSubdivisionCandidate &splitCandidate)
{
	VspBspTraversalData frontData;
	VspBspTraversalData backData;
	
	BspLeaf *leaf = dynamic_cast<BspLeaf *>(splitCandidate.mParentData.mNode);
	
	// compute locally best split plane
    const bool costRatioViolated = 
		SelectPlane(splitCandidate.mSplitPlane, 
					leaf, 
					splitCandidate.mParentData, 
					frontData, 
					backData, 
					splitCandidate.mSplitAxis);

	// max cost threshold violated?
	splitCandidate.mMaxCostMisses = costRatioViolated ? 
		splitCandidate.mParentData.mMaxCostMisses : 
		splitCandidate.mParentData.mMaxCostMisses + 1;

	float oldRenderCost;

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

	splitCandidate.mRenderCostDecr = renderCostDecr;

	// TODO: geometry could be reused
	delete frontData.mGeometry;
	delete backData.mGeometry;

	// set priority for queue
#if 0
	const float priority = (float)-data.mDepth;
#else	

	// take render cost of node into account 
	// otherwise danger of being stuck in a local minimum!!
	const float factor = mRenderCostDecreaseWeight;
	const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;
#endif
	
	splitCandidate.mPriority = priority;
}


void VspBspTree::EvalSubdivisionStats(const VspBspTraversalData &tData,
									  const VspBspTraversalData &tFrontData,
									  const VspBspTraversalData &tBackData)
{
	const float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
	const float cBack = (float)tBackData.mPvs * tBackData.mProbability;
	const float cData = (float)tData.mPvs * tData.mProbability;
	
	const float costDecr = 
		(cFront + cBack - cData) / mBoundingBox.GetVolume();

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

	AddSubdivisionStats(mBspStats.Leaves(),
						-costDecr,
						0,
						mTotalCost,
						(float)mTotalPvsSize / (float)mBspStats.Leaves());
}


BspInterior *VspBspTree::SubdivideNode(const Plane3 &splitPlane,
									   VspBspTraversalData &tData,
									   VspBspTraversalData &frontData,
									   VspBspTraversalData &backData,
									   PolygonContainer &coincident)
{
	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
	
	//-- 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(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
	
	// 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,
									   splitPlane,
									   mBoundingBox,
									   //0.0f);
									   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;

			// should never come here: wrong volume !!!
			if (0)
			{
				if (frontData.mProbability < -0.00001)
					Debug << "fatal error f: " << frontData.mProbability << endl;
				if (backData.mProbability < -0.00001)
					Debug << "fatal error b: " << backData.mProbability << endl;

				// clamp because of precision issues
				if (frontData.mProbability < 0)	frontData.mProbability = 0;
				if (backData.mProbability < 0) backData.mProbability = 0;
			}
		}
	}

	
    // subdivide polygons
	SplitPolygons(splitPlane,
				  *tData.mPolygons,
	              *frontData.mPolygons,
				  *backData.mPolygons,
				  coincident);



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

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

	mBspStats.nodes += 2;

    
	BspInterior *interior = new BspInterior(splitPlane);

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


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

	interior->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();
  
	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;
		}
	  
		// only count termination objects?
		if (COUNT_ORIGIN_OBJECTS && ray->mOriginObject) 
		{
			if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
				madeContrib = true;

			sc += contribution;
		}
		
		sampleContributions += sc;
		
		if (madeContrib)
			++ contributingSamples;
	}
}


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

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

	mLocalSubdivisionCandidates->reserve(requestedSize);

	if (0)
	{	// float values => don't compare with exact values
		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);
		float pos = (*ri).ExtrapOrigin(axis);

		// clamp to min / max band
		if (0) ClipValue(pos, minBand, maxBand);
		
		mLocalSubdivisionCandidates->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);

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

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


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

	if (mMaxTests < (int)rays.size())
	{
		GetRayInfoSets(rays, mMaxTests, usedRays);
	}
	else
	{
		usedRays = rays;
	}

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

	const float sizeBox = maxBox - minBox;

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

	SortSubdivisionCandidates(usedRays, axis, minBand, maxBand);

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

	float pvsl = 0;
	float pvsr = (float)pvsSize;

	float pvsBack = pvsl;
	float 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 = 99999999.0f;
	bool splitPlaneFound = false;

	Intersectable::NewMail();
   	RayInfoContainer::const_iterator ri, ri_end = usedRays.end();

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

		if (COUNT_ORIGIN_OBJECTS && 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 = mLocalSubdivisionCandidates->end();

	for (ci = mLocalSubdivisionCandidates->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 (COUNT_ORIGIN_OBJECTS && oObject && !oObject->Mailed())
					{
						oObject->Mail();
						++ pvsl;
					}

					if (tObject && !tObject->Mailed())
					{
						tObject->Mail();
						++ pvsl;
					}

					break;
				}
			case SortableEntry::ERayMax:
				{
					if (COUNT_ORIGIN_OBJECTS && oObject)
					{
						if (-- oObject->mCounter == 0)
							-- pvsr;
					}

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

					break;
				}
		}
		
		
		// Note: we compare size of bounding boxes of front and back side because
		// of efficiency reasons (otherwise a new geometry would have to be computed
		// in each step and incremential evaluation would be difficult.
		// but then errors happen if the geometry is not an axis aligned box 
		// (i.e., if a geometry aligned split was taken before)
		// question: is it sufficient to make this approximation?
		if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
		{
			sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);

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

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

			if (sum < minSum)
			{
				splitPlaneFound = true;

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

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

	const float pOverall = sizeBox;

	const float pBack = position - minBox;
	const float pFront = maxBox - position;
	
	const float penaltyOld = EvalPvsPenalty((float)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 = mPvsFactor * 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 VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
										 const VspBspTraversalData &tData,
										 int &axis,
										 BspNodeGeometry **frontGeom,
										 BspNodeGeometry **backGeom,
										 float &pFront,
										 float &pBack,
										 const bool isKdNode)
{
	float nPosition[3];
	float nCostRatio[3];
	float nProbFront[3];
	float nProbBack[3];

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

	// set to NULL, so I can find out which gemetry was stored
	for (int i = 0; i < 3; ++ i)
	{
		nFrontGeom[i] = NULL;
		nBackGeom[i] = NULL;
	}

	// create bounding box of node geometry
	AxisAlignedBox3 box;
		
	//TODO: for kd split geometry already is box => only take minmax vertices
	if (1)
	{	// get bounding box from geometry
		tData.mGeometry->GetBoundingBox(box);
	}
	else
	{
		box.Initialize();
		RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();

		for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
			box.Include((*ri).ExtrapTermination());
	}


	int sAxis = 0;
	int bestAxis;

	// if max cost ratio is exceeded, take split along longest axis instead
	const float maxCostRatioForArbitraryAxis = 0.9f;

	if (mUseDrivingAxisIfMaxCostViolated)
		bestAxis = box.Size().DrivingAxis();
	else
		bestAxis = -1;

#if 0
	// maximum cost ratio for axis to be valid:
	// if exceeded, spatial mid split is used instead
	const maxCostRatioForHeur = 0.99f;
#endif

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

		
	//Debug << "use special axis: " << useSpecialAxis << endl;
	//Debug << "axis: " << sAxis << " drivingaxis: " << box.Size().DrivingAxis();
	
	for (axis = 0; axis < 3 ; ++ axis)
	{
		if (!useSpecialAxis || (axis == sAxis))
		{
			if (mUseCostHeuristics)
			{
				//-- place split plane using heuristics
				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;
				Vector3 normal(0,0,0); normal[axis] = 1.0f;
				
				// allows faster split because we have axis aligned kd tree boxes
				if (isKdNode)
				{
					nCostRatio[axis] = EvalAxisAlignedSplitCost(tData,
																box,
																axis,
																nPosition[axis],
																nProbFront[axis], 
																nProbBack[axis]);
					
					// create back geometry from box

					// NOTE: the geometry is returned from the function so we
					// don't have to recompute it when possible
					Vector3 pos;
					
					pos = box.Max(); pos[axis] = nPosition[axis];
					AxisAlignedBox3 bBox(box.Min(), pos);
					
					PolygonContainer fPolys;
					bBox.ExtractPolys(fPolys);

					nBackGeom[axis] = new BspNodeGeometry(fPolys);
	
					////////////
					//-- create front geometry from box

					pos = box.Min(); pos[axis] = nPosition[axis];
					AxisAlignedBox3 fBox(pos, box.Max());

					PolygonContainer bPolys;
					fBox.ExtractPolys(bPolys);
					nFrontGeom[axis] = new BspNodeGeometry(bPolys);
				}
				else
				{
					nFrontGeom[axis] = new BspNodeGeometry();
					nBackGeom[axis] = new BspNodeGeometry();

					nCostRatio[axis] = 
						EvalSplitPlaneCost(Plane3(normal, nPosition[axis]), 
										   tData, *nFrontGeom[axis], *nBackGeom[axis],
										   nProbFront[axis], nProbBack[axis]);
				}
			}
						
			
			if (mUseDrivingAxisIfMaxCostViolated)
			{
				// we take longest axis split if cost ratio exceeds threshold
				if (nCostRatio[axis] < min(maxCostRatioForArbitraryAxis, nCostRatio[bestAxis]))
				{
					bestAxis = axis;
				}
				/*else if (nCostRatio[axis] < nCostRatio[bestAxis])
				{
					Debug << "taking split along longest axis (" << bestAxis << ") instead of  (" << axis << ")" << endl;
				}*/

			}
			else
			{
				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]);

	//Debug << "best axis: " << bestAxis << " pos " << nPosition[bestAxis] << endl;

	return nCostRatio[bestAxis];
}


bool VspBspTree::SelectPlane(Plane3 &bestPlane,
							 BspLeaf *leaf,
							 VspBspTraversalData &data,							 
							 VspBspTraversalData &frontData,
							 VspBspTraversalData &backData,
							 int &splitAxis)
{
	// HACK matt: subdivide regularily to certain depth
	if (data.mDepth < 0)  	// question matt: why depth < 0 ?
	{
		cout << "depth: " << data.mDepth << endl;

		// return axis aligned split
		AxisAlignedBox3 box;
		box.Initialize();
	
		// create bounding box of region
		data.mGeometry->GetBoundingBox(box);
	
		const int axis = box.Size().DrivingAxis();
		const Vector3 position = (box.Min()[axis] + box.Max()[axis]) * 0.5f;

		Vector3 norm(0,0,0); norm[axis] = 1.0f;
		bestPlane = Plane3(norm, position);
		splitAxis = axis;

		return true;
	}

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


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

	candidateCost = 99999999.0f;

	// as a variant, we take axis aligned split only if there is 
	// more polygon available to guide the split
	if (!mUsePolygonSplitIfAvailable || data.mPolygons->empty())
	{
		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);
	}
    
#ifdef GTP_DEBUG
	Debug << "plane lowest cost: " << lowestCost << endl;
#endif

	// exeeded relative max cost ratio
	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::EvalRenderCostDecrease(const Plane3 &candidatePlane,
										 const VspBspTraversalData &data,
										 float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;
	float totalPvs = 0;

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


	// create unique ids for pvs heuristics
	GenerateUniqueIdsForPvs();
	
	for (int i = 0; i < data.mRays->size(); ++ i)
	{
		RayInfo rayInf = (*data.mRays)[i];

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

		// find front and back pvs for origing and termination object
		AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);

		if (COUNT_ORIGIN_OBJECTS)
		{
			AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
		}
	}


	BspNodeGeometry geomFront;
	BspNodeGeometry geomBack;

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

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

		// something is wrong with the volume
		if (0 && ((pFront < 0.0) || (pBack < 0.0)))
		{
			Debug << "ERROR in volume:\n"
				  << "volume f :" << pFront << " b: " << pBack << " p: " << pOverall 
				  << ", real volume f: " << pFront << " b: " << geomBack.GetVolume()
				  << ", #polygons f: " << geomFront.Size() << " b: " << geomBack.Size() << " p: " << data.mGeometry->Size() << endl;
		}
	}
	else
	{
		pFront = geomFront.GetArea();
		pBack = geomBack.GetArea();
	}
	

	// -- pvs rendering heuristics

	// upper and lower bounds
	const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
	const float upperPvsLimit = (float)mViewCellsManager->GetMaxPvsSize();

	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;

	const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBoundingBox.GetVolume();
	
	// take render cost of node into account to avoid being stuck in a local minimum
	normalizedOldRenderCost = oldRenderCost / mBoundingBox.GetVolume();
	
	return renderCostDecrease;
}


float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
									 const VspBspTraversalData &data,
									 BspNodeGeometry &geomFront,
									 BspNodeGeometry &geomBack,
									 float &pFront,
									 float &pBack) const
{
	float totalPvs = 0;
	float pvsFront = 0;
	float pvsBack = 0;
	
	// overall probability is used as normalizer
	float pOverall = 0;

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

	int numTests; // the number of tests

	// if random samples shold be taken instead of testing all the rays
	bool useRand; 

	if ((int)data.mRays->size() > mMaxTests)
	{
		useRand = true;
		numTests = mMaxTests;
	}
	else
	{
		useRand = false;
		numTests = (int)data.mRays->size();
	}
	
	// create unique ids for pvs heuristics
	GenerateUniqueIdsForPvs();

	for (int i = 0; i < numTests; ++ 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);

		// find front and back pvs for origing and termination object
		AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);

		if (COUNT_ORIGIN_OBJECTS)
		{
			AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
		}
	}

	// construct child geometry with regard to the candidate split plane
	bool splitSuccessFull = data.mGeometry->SplitGeometry(geomFront,
														  geomBack,
														  candidatePlane,
														  mBoundingBox,
														  //0.0f);
														  mEpsilon);

	pOverall = data.mProbability;

	if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
	{
		pFront = geomFront.GetVolume();
		pBack = pOverall - pFront;
		
		// HACK: precision issues possible for unbalanced split => don't take this split!
		if (1 && 
			(!splitSuccessFull || (pFront <= 0) || (pBack <= 0) || 
			!geomFront.Valid() || !geomBack.Valid()))
		{
			//Debug << "error f: " << pFront << " b: " << pBack << endl;

			// high penalty for degenerated / wrong split
			return 99999.9f;
		}
	}
	else
	{
		pFront = geomFront.GetArea();
		pBack = geomBack.GetArea();
	}
	
	////////
	//-- pvs rendering heuristics

	const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
	const float upperPvsLimit = (float)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;
	}

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

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

	return cost;
}


int VspBspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box, 
										ViewCellContainer &viewCells) const
{
	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();

		const int side = geom->ComputeIntersection(box);
		
		switch (side)
		{
		case -1:
			// node geometry is contained in box
			CollectViewCells(node, true, viewCells, true);
			break;

		case 0:
			if (node->IsLeaf())
			{
				BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
			
				if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
				{
					leaf->GetViewCell()->Mail();
					viewCells.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(),
									mBoundingBox,
									//0.0000001f);
									mEpsilon);

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

	return (int)viewCells.size();
}


float VspBspTree::EvalAxisAlignedSplitCost(const VspBspTraversalData &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);

		if (COUNT_ORIGIN_OBJECTS)
		{
			AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
		}
	}


	//-- pvs heuristics

	float pOverall = data.mProbability;

	// note: we use a simplified computation assuming that we always do a
	// spatial mid split	
	
	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 GTP_DEBUG
	Debug << "axis: " << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
	Debug << "p: " << pFront << " " << pBack << " " << pOverall << endl;
#endif

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

	return  (mCtDivCi + newCost) / oldCost;
}


inline void VspBspTree::AddObjToPvs(Intersectable *obj,
										 const int cf,
										 float &frontPvs,
										 float &backPvs,
										 float &totalPvs) const
{
	if (!obj)
		return;
#if 0
	const float renderCost = mViewCellsManager->EvalRenderCost(obj);
#else
	const int renderCost = 1;
#endif
	// 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 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().EvalPvsCost() <= 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 mBoundingBox;
}


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


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

	mBspStats.pvs += data.mPvs;

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

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


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

	stack<BspRayTraversalData> tQueue;

	float maxt, mint;

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

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

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

	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
			{
				// matt: 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::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())
			{
				ViewCellLeaf *leafVc = dynamic_cast<BspLeaf *>(node)->GetViewCell();

				ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
						
				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());
		}
	}
}


void VspBspTree::CollapseViewCells()
{
// TODO
#if HAS_TO_BE_REDONE
	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::PreprocessPolygons(PolygonContainer &polys)
{
	// preprocess: throw out polygons coincident to the view space box (not needed)
	PolygonContainer boxPolys;
	
	mBoundingBox.ExtractPolys(boxPolys);
	vector<Plane3> boxPlanes;

	PolygonContainer::iterator pit, pit_end = boxPolys.end();

	// extract planes of box
	// TODO: can be done more elegantly than first extracting polygons
	// and take their planes
	for (pit = boxPolys.begin(); pit != pit_end; ++ pit)
	{
		boxPlanes.push_back((*pit)->GetSupportingPlane());
	}

	pit_end = polys.end();

	for (pit = polys.begin(); pit != pit_end; ++ pit)
	{
		vector<Plane3>::const_iterator bit, bit_end = boxPlanes.end();
		
		for (bit = boxPlanes.begin(); (bit != bit_end) && (*pit); ++ bit)
		{
			const int cf = (*pit)->ClassifyPlane(*bit, mEpsilon);

			if (cf == Polygon3::COINCIDENT)
			{
				DEL_PTR(*pit);
				//Debug << "coincident!!" << endl;
			}
		}
	}

	// remove deleted entries after swapping them to end of vector
	for (int i = 0; i < (int)polys.size(); ++ i)
	{
		while (!polys[i] && (i < (int)polys.size()))
		{
			swap(polys[i], polys.back());
			polys.pop_back();
		}
	}
}


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

				++ splits;

				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->GetBack() != 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;
	vector<Plane3> candidatePlanes;

	vector<Plane3>::const_iterator pit, pit_end = halfSpaces.end();

	// bounded planes are added to the polygons
	for (pit = halfSpaces.begin(); pit != pit_end; ++ pit)
	{
		Polygon3 *p = GetBoundingBox().CrossSection(*pit);

		if (p->Valid(mEpsilon))
		{
			candidatePolys.push_back(p);
			candidatePlanes.push_back(*pit);
		}
	}

	// 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(mBoundingBox.GetFace(i).mVertices[j]);
		}

		Polygon3 *poly = new Polygon3(vertices);

		candidatePolys.push_back(poly);
		candidatePlanes.push_back(poly->GetSupportingPlane());
	}

	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 (backPoly->Valid(mEpsilon))
					{
						candidatePolys[i] = backPoly;
					}
					else
					{
						DEL_PTR(backPoly);
					}

					// outside, don't need this
					DEL_PTR(frontPoly);
					break;

				// polygon outside of halfspace
				case Polygon3::FRONT_SIDE:
					DEL_PTR(candidatePolys[i]);
					break;

				// just take polygon as it is
				case Polygon3::BACK_SIDE:
				case Polygon3::COINCIDENT:
				default:
					break;
			}
		}

		if (candidatePolys[i])
		{
			geom.Add(candidatePolys[i], candidatePlanes[i]);
		}
	}
}


bool VspBspTree::IsOutOfBounds(ViewCell *vc) const
{
	return vc->GetId() == OUT_OF_BOUNDS_ID;
}


void VspBspTree::SetViewCellsTree(ViewCellsTree *viewCellsTree)
{
	mViewCellsTree = viewCellsTree;
}


void VspBspTree::ConstructGeometry(ViewCell *vc, 
								   BspNodeGeometry &vcGeom) const
{
	// if false, cannot construct geometry for interior leaf
	if (!mViewCellsTree)
		return;

	ViewCellContainer leaves;
	mViewCellsTree->CollectLeaves(vc, leaves);

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

	for (it = leaves.begin(); it != it_end; ++ it)
	{
		if (IsOutOfBounds(*it))
			continue;
		
		BspViewCell *bspVc = dynamic_cast<BspViewCell *>(*it);
		vector<BspLeaf *>::const_iterator bit, bit_end = bspVc->mLeaves.end();

		for (bit = bspVc->mLeaves.begin(); bit != bit_end; ++ bit)
		{
			BspLeaf *l = *bit;
			ConstructGeometry(l, vcGeom);
		}
	}
}


int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
							  const bool onlyUnmailed) const
{
	stack<bspNodePair> nodeStack;
	
	BspNodeGeometry nodeGeom;
	ConstructGeometry(n, nodeGeom);
//	const float eps = 0.5f;
	const float eps = 0.01f;
	// 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->GetPolys(),
													halfSpaces[i],
													eps);

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

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

						if (cf == Polygon3::FRONT_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.GetPolys(),
												   interior->GetPlane(),
												   eps);
			
			BspNode *front = interior->GetFront();
			BspNode *back = interior->GetBack();
            
			BspNodeGeometry *fGeom = new BspNodeGeometry();
			BspNodeGeometry *bGeom = new BspNodeGeometry();

			geom->SplitGeometry(*fGeom,
								*bGeom,
								interior->GetPlane(),
								mBoundingBox,
								//0.0000001f);
								eps);
		
			if (cf == Polygon3::BACK_SIDE)
			{
				nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
				DEL_PTR(fGeom);
			}
			else
			{
				if (cf == Polygon3::FRONT_SIDE)
				{
					nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
					DEL_PTR(bGeom);
				}
				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);
	
	float eps = 0.01f;
	// 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.GetPolys(),
												   interior->GetPlane(),
												   eps);
			
			BspNode *front = interior->GetFront();
			BspNode *back = interior->GetBack();
            
			BspNodeGeometry *fGeom = new BspNodeGeometry();
			BspNodeGeometry *bGeom = new BspNodeGeometry();

			geom->SplitGeometry(*fGeom,
								*bGeom,
								interior->GetPlane(),
								mBoundingBox,
								//0.0000001f);
								eps);
		
			if (cf == Polygon3::BACK_SIDE)
			{
				nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
				DEL_PTR(fGeom);
				}
			else
			{
				if (cf == Polygon3::FRONT_SIDE)
				{
					nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
					DEL_PTR(bGeom);
				}
				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.GetPolys(), 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 (COUNT_ORIGIN_OBJECTS && 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,
								ViewCellContainer &viewcells)
{
	int hits = 0;
	stack<BspRayTraversalData> tStack;

	float mint = 0.0f, maxt = 1.0f;

	//ViewCell::NewMail();

	Vector3 entp = origin;
	Vector3 extp = termination;

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

	float t;
	const float thresh = 1e-6f; // matt: change this to adjustable value
	
	while (1)
	{
		if (!node->IsLeaf())
		{
			BspInterior *in = dynamic_cast<BspInterior *>(node);

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

			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 // one of the ray end points is on the plane
			{	
				// NOTE: what to do if ray is coincident with plane?
				if (extSide < 0)
					node = in->GetBack();
				else //if (extSide > 0)
					node = in->GetFront();
				//else break; // coincident => count no intersections

				continue; // no far child
			}

			// push data for far child
			tStack.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;
			
			// question: always contribute to leaf or to currently active view cell?
			if (0)
				viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
			else
				viewCell = leaf->GetViewCell();

			if (!viewCell->Mailed())
			{
				viewcells.push_back(viewCell);
				viewCell->Mail();
				++ hits;
			}

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

			entp = extp;
			
			const BspRayTraversalData &s = tStack.top();

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

			tStack.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 HAS_TO_BE_REDONE
	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 HAS_TO_BE_REDONE
	(void)CollapseTree(mRoot, collapsed);

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


void VspBspTree::RepairViewCellsLeafLists()
{
// TODO
#if HAS_TO_BE_REDONE
	// 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
}


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;
		geom->GetBoundingBox(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(),
									mBoundingBox,
									//0.0000001f);
									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;
		
		// appoximate neighbor search has slightl relaxed constraints
		if (1)
			FindNeighbors(leaf, neighbors, true);
		else
			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());

				if (!leaf->GetViewCell()->GetPvs().Empty() || 
					!(*nit)->GetViewCell()->GetPvs().Empty() ||
                    leaf->IsSibling(*nit))
				{
					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;

		iit = ray->mViewCells.begin();
		BspViewCell *bspVc = dynamic_cast<BspViewCell *>(*(iit ++));
		BspLeaf *leaf = bspVc->mLeaves[0];
		
		// 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->mLeaves[0]; // exactly one leaf

			// 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, const bool active)
{
	if (mRoot == NULL)
		return NULL;

	stack<BspNode *> nodeStack;
	nodeStack.push(mRoot);
  
	ViewCellLeaf *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());
		}
	}
  
	if (active)
		return mViewCellsTree->GetActiveViewCell(viewcell);
	else
		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(OUT_STREAM &stream)
{
	ExportNode(mRoot, stream);
	return true;
}


void VspBspTree::ExportNode(BspNode *node, OUT_STREAM &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;
	}
}

}