#include "Plane3.h"
#include "ViewCellBsp.h"
#include "Mesh.h"
#include "common.h"
#include "ViewCell.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include "Triangle3.h"
#include "Tetrahedron3.h"
#include "ViewCellsManager.h"
#include "Exporter.h"
#include "Plane3.h"
#include <stack>

namespace GtpVisibilityPreprocessor {

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

int BspNode::sMailId = 1;

/** Evaluates split plane classification with respect to the plane's 
	contribution for a minimum number splits in the tree.
*/
const float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
/** Evaluates split plane classification with respect to the plane's 
	contribution for a balanced tree.
*/
const float BspTree::sBalancedPolysTable[] = {1, -1, 0, 0};

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

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


/******************************************************************/
/*                  class BspNode implementation                  */
/******************************************************************/


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


BspNode::BspNode(BspInterior *parent): 
mParent(parent), mTreeValid(true)
{}


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


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

	return depth;
}


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


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


/****************************************************************/
/*              class BspInterior implementation                */
/****************************************************************/


BspInterior::BspInterior(const Plane3 &plane): 
mPlane(plane), mFront(NULL), mBack(NULL)
{}

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


void BspInterior::SetupChildLinks(BspNode *b, BspNode *f) 
{
    mBack = b;
    mFront = f;
}


/****************************************************************/
/*                  class BspLeaf implementation                */
/****************************************************************/


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


BspLeaf::~BspLeaf()
{
	DEL_PTR(mPvs);

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

		if (!ray->IsActive())
			delete ray;
	}
}


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


BspLeaf::BspLeaf(BspInterior *parent): 
BspNode(parent), mViewCell(NULL), mPvs(NULL)
{}



BspLeaf::BspLeaf(BspInterior *parent, ViewCellLeaf *viewCell): 
BspNode(parent), mViewCell(viewCell), mPvs(NULL)
{
}


/*********************************************************************/
/*                       class BspTree implementation                */
/*********************************************************************/

BspTree::BspTree():  
mRoot(NULL), 
mUseAreaForPvs(false),
mUsePredefinedViewCells(false),
mTimeStamp(1),
mViewCellsTree(NULL),
mOutOfBoundsCellPartOfTree(false),
mOutOfBoundsCell(NULL)
{
	Randomize(); // initialise random generator for heuristics
	mOutOfBoundsCell = GetOrCreateOutOfBoundsCell();

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

	Environment::GetSingleton()->GetIntValue("BspTree.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.minPvs", mTermMinPvs);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.minPolygons", mTermMinPolys);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.minRays", mTermMinRays);
	Environment::GetSingleton()->GetFloatValue("BspTree.Termination.minProbability", mTermMinProbability);	
	Environment::GetSingleton()->GetFloatValue("BspTree.Termination.maxRayContribution", mTermMaxRayContribution);
	Environment::GetSingleton()->GetFloatValue("BspTree.Termination.minAccRayLenght", mTermMinAccRayLength);

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

	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.verticalSplits", mVerticalSplitsFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.largestPolyArea", mLargestPolyAreaFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.blockedRays", mBlockedRaysFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.leastRaySplits", mLeastRaySplitsFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.balancedRays", mBalancedRaysFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.pvs", mPvsFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.leastSplits" , mLeastSplitsFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.balancedPolys", mBalancedPolysFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Factor.balancedViewCells", mBalancedViewCellsFactor);
	Environment::GetSingleton()->GetFloatValue("BspTree.Termination.ct_div_ci", mCtDivCi);


	///////////
	//-- termination criteria for axis aligned split

	Environment::GetSingleton()->GetFloatValue("BspTree.Termination.AxisAligned.ct_div_ci", mAxisAlignedCtDivCi);
	Environment::GetSingleton()->GetFloatValue("BspTree.Termination.maxCostRatio", mMaxCostRatio);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.AxisAligned.minPolys", 
							 mTermMinPolysForAxisAligned);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.AxisAligned.minRays", 
							 mTermMinRaysForAxisAligned);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.AxisAligned.minObjects", 
							 mTermMinObjectsForAxisAligned);
	
	
	//////////////
	//-- partition criteria

	Environment::GetSingleton()->GetIntValue("BspTree.maxPolyCandidates", mMaxPolyCandidates);
	Environment::GetSingleton()->GetIntValue("BspTree.maxRayCandidates", mMaxRayCandidates);
	Environment::GetSingleton()->GetIntValue("BspTree.splitPlaneStrategy", mSplitPlaneStrategy);
	Environment::GetSingleton()->GetFloatValue("BspTree.AxisAligned.splitBorder", mSplitBorder);
	Environment::GetSingleton()->GetIntValue("BspTree.maxTests", mMaxTests);
	Environment::GetSingleton()->GetIntValue("BspTree.Termination.maxViewCells", mMaxViewCells);

	Environment::GetSingleton()->GetFloatValue("BspTree.Construction.epsilon", mEpsilon);

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

	/////////////////////7

	Debug << "BSP options: " << endl;
    Debug << "BSP max depth: " << mTermMaxDepth << endl;
	Debug << "BSP min PVS: " << mTermMinPvs << endl;
	Debug << "BSP min probability: " << mTermMinProbability << endl;
	Debug << "BSP max polys: " << mTermMinPolys << endl;
	Debug << "BSP max rays: " << mTermMinRays << endl;
	Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
	Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;

	Debug << "Split plane strategy: ";
	if (mSplitPlaneStrategy & RANDOM_POLYGON)
		Debug << "random polygon ";
	if (mSplitPlaneStrategy & AXIS_ALIGNED)
		Debug << "axis aligned ";
	if (mSplitPlaneStrategy & LEAST_SPLITS)     
		Debug << "least splits ";
	if (mSplitPlaneStrategy & BALANCED_POLYS)
		Debug << "balanced polygons ";
	if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
		Debug << "balanced view cells ";
	if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
		Debug << "largest polygon area ";
	if (mSplitPlaneStrategy & VERTICAL_AXIS)
		Debug << "vertical axis ";
	if (mSplitPlaneStrategy & BLOCKED_RAYS)
		Debug << "blocked rays ";
	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
		Debug << "least ray splits ";
	if (mSplitPlaneStrategy & BALANCED_RAYS)
		Debug << "balanced rays ";
	if (mSplitPlaneStrategy & PVS)
		Debug << "pvs";

	Debug << endl;
}


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


const BspTreeStatistics &BspTree::GetStatistics() const 
{
	return mStat;
}


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

	return mOutOfBoundsCell;
}

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

#ifdef _Debug
	Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
#endif

	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:
				{
					Polygon3 *front_piece = new Polygon3(poly->mParent);
					Polygon3 *back_piece = new Polygon3(poly->mParent);

					//-- split polygon into front and back part
					poly->Split(plane, 
								*front_piece, 
								*back_piece,
								mEpsilon);
					
					++ splits; // increase number of splits

					//-- inherit rays from parent polygon for blocked ray criterium
					poly->InheritRays(*front_piece, *back_piece);
				
					// check if polygons still valid 
					if (front_piece->Valid(mEpsilon))
						frontPolys.push_back(front_piece);
					else
						DEL_PTR(front_piece);
				
					if (back_piece->Valid(mEpsilon))
						backPolys.push_back(back_piece);
					else				
						DEL_PTR(back_piece);
				
#ifdef GTP_DEBUG
					Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
#endif
					DEL_PTR(poly);
				}
				break;
			default:
                Debug << "SHOULD NEVER COME HERE\n";
				break;
		}
	}

	return splits;
}


void BspTreeStatistics::Print(ostream &app) const
{
	app << "===== BspTree statistics ===============\n";

	app << setprecision(4);

	app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";

	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";

	app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";

	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";

	app << "#N_POLYSPLITS ( Number of polygon splits )\n" << polySplits << "\n";

	app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits[0] + splits[1] + splits[2] << endl;

	app << "#N_SPLITS ( Number of splits in axes x y z)\n";

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

	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
		<<	maxDepthNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n" 
		<< minPvsNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n" 
		<< minRaysNodes * 100 / (double)Leaves() << endl;

	app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
		<< maxCostNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
		<< minProbabilityNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n" 
		<<	maxRayContribNodes * 100 / (double)Leaves() << endl;

	app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;

	app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;

	app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;

	app << "#N_INPUTPOLYGONS (number of input polygons )\n" << polys << endl;

	app << "#N_INVALIDLEAVES (number of invalid leaves )\n" << invalidLeaves << endl;

	app << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
	//app << "#N_PVS: " << pvs << endl;

	app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
		 (polys + polySplits) / (double)polys << endl;
	
	app << "===== END OF BspTree statistics ==========\n";
}


BspTree::~BspTree()
{
	DEL_PTR(mRoot);
	
	if (mOutOfBoundsCellPartOfTree)
	{
		// out of bounds cell not part of tree =>
		// delete manually
		DEL_PTR(mOutOfBoundsCell);
	}
}


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


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


void BspTree::InsertViewCell(ViewCellLeaf *viewCell)
{
	// don't generate new view cell, insert this view cell
	mUsePredefinedViewCells = true;
	PolygonContainer *polys = new PolygonContainer();

	// extract polygons that guide the split process
	mStat.polys += AddMeshToPolygons(viewCell->GetMesh(), *polys, viewCell);
	mBbox.Include(viewCell->GetBox()); // add to BSP aabb

	InsertPolygons(polys);
}


void BspTree::InsertPolygons(PolygonContainer *polys)
{	
	BspTraversalStack tStack;

	// traverse existing tree or create new tree
    if (!mRoot)
		mRoot = new BspLeaf();

	tStack.push(BspTraversalData(mRoot, 
								 polys, 
								 0, 
								 mOutOfBoundsCell, 
								 new BoundedRayContainer(), 
								 0, 
								 mUseAreaForPvs ? mBbox.SurfaceArea() : mBbox.GetVolume(), 
								 new BspNodeGeometry()));

	while (!tStack.empty())
	{
		// filter polygons donw the tree
		BspTraversalData tData = tStack.top();
	    tStack.pop();
			
		if (!tData.mNode->IsLeaf())
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(tData.mNode);

			///////////////////
			//-- filter view cell polygons down the tree until a leaf is reached
			if (!tData.mPolygons->empty())
			{
				PolygonContainer *frontPolys = new PolygonContainer();
				PolygonContainer *backPolys = new PolygonContainer();
				PolygonContainer coincident;

				int splits = 0;
		
				// split viewcell polygons with respect to split plane
				splits += SplitPolygons(interior->GetPlane(),
										*tData.mPolygons, 
										*frontPolys, 
										*backPolys, 
										coincident);
				
				// extract view cells associated with the split polygons
				ViewCellLeaf *frontViewCell = mOutOfBoundsCell;
				ViewCellLeaf *backViewCell = mOutOfBoundsCell;
			
				BspTraversalData frontData(interior->GetFront(), 
										   frontPolys, 
										   tData.mDepth + 1,
										   mOutOfBoundsCell,	
										   tData.mRays,
										   tData.mPvs,
										   mUseAreaForPvs ? mBbox.SurfaceArea() : mBbox.GetVolume(), 
										   new BspNodeGeometry());

				BspTraversalData backData(interior->GetBack(), 
										  backPolys,
										  tData.mDepth + 1,
										  mOutOfBoundsCell,	
										  tData.mRays,
										  tData.mPvs,
										   mUseAreaForPvs ? mBbox.SurfaceArea() : mBbox.GetVolume(), 
										  new BspNodeGeometry());

				if (mUsePredefinedViewCells)
				{
					ExtractViewCells(frontData,
									 backData,
									 coincident,
									 interior->mPlane);
				}

				// don't need coincident polygons anymore
				CLEAR_CONTAINER(coincident);

				mStat.polySplits += splits;

				// push the children on the stack
				tStack.push(frontData);
				tStack.push(backData);
			}

			// cleanup
			DEL_PTR(tData.mPolygons);
			DEL_PTR(tData.mRays);
		}
		else
		{
			// reached leaf => subdivide current viewcell
			BspNode *subRoot = Subdivide(tStack, tData);
		}
	}
}


int BspTree::AddMeshToPolygons(Mesh *mesh, 
							   PolygonContainer &polys, 
							   MeshInstance *parent)
{
	FaceContainer::const_iterator fi, fi_end = mesh->mFaces.end();
	
	// copy the face data to polygons
	for (fi = mesh->mFaces.begin(); fi != fi_end; ++ fi)
	{
		Polygon3 *poly = new Polygon3((*fi), mesh);
		
		if (poly->Valid(mEpsilon))
		{
			poly->mParent = parent; // set parent intersectable
			polys.push_back(poly);
		}
		else
		{
			DEL_PTR(poly);
		}
	}
	return (int)mesh->mFaces.size();
}


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

	for (int i = 0; i < limit; ++ i)
	{
		Mesh *mesh = viewCells[i]->GetMesh();

		if (mesh) 
		{
			// copy the mesh into polygons and add to BSP tree aabb
			mBbox.Include(viewCells[i]->GetBox()); 
			polysSize += AddMeshToPolygons(mesh, polys, viewCells[i]);
		}
	}

	return polysSize;
}


int BspTree::AddToPolygonSoup(const ObjectContainer &objects,
							  PolygonContainer &polys, 
							  int maxObjects,
							  bool addToBbox)
{
	int limit = (maxObjects > 0) ? 
		Min((int)objects.size(), maxObjects) : (int)objects.size();
  
	for (int i = 0; i < limit; ++i)
	{
		Intersectable *object = objects[i];
		Mesh *mesh = NULL;

		switch (object->Type()) // extract the meshes
		{
		case Intersectable::MESH_INSTANCE:
			mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
			break;
		case Intersectable::VIEW_CELL:
			mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
			break;
		case Intersectable::TRANSFORMED_MESH_INSTANCE:
			{
				TransformedMeshInstance *mi = dynamic_cast<TransformedMeshInstance *>(object);

				if (!mi->GetMesh())	
					break;
				mesh = new Mesh();
				mi->GetTransformedMesh(*mesh);
				break;
			}
		default:
			Debug << "intersectable type not supported" << endl;
			break;
		}
		
        if (!mesh) continue;

		// copy the mesh data to polygons
		if (addToBbox)
		{
			mBbox.Include(object->GetBox()); // add to BSP tree aabb
		}

		AddMeshToPolygons(mesh, polys, mOutOfBoundsCell);

		// cleanup
		if (object->Type() == Intersectable::TRANSFORMED_MESH_INSTANCE)
		{
			DEL_PTR(mesh);
		}	
	}

	return (int)polys.size();
}

	
void BspTree::Construct(const ViewCellContainer &viewCells)
{
	// construct hierarchy over the given view cells
	mUsePredefinedViewCells = true;

	mStat.nodes = 1;
	mBbox.Initialize(); // initialise bsp tree bounding box

 	// copy view cell meshes into one big polygon soup
	PolygonContainer *polys = new PolygonContainer();
	mStat.polys = AddToPolygonSoup(viewCells, *polys);

	Exporter *expo = Exporter::GetExporter("dummy2.wrl");
	expo->ExportPolygons(*polys);
	delete expo;
	// construct tree from the view cell polygons
	Construct(polys, new BoundedRayContainer());
}


void BspTree::Construct(const ObjectContainer &objects)
{
	// generate new view cells for this type
	mUsePredefinedViewCells = false;

	mStat.nodes = 1;
	mBbox.Initialize(); // initialise bsp tree bounding box
	
	PolygonContainer *polys = new PolygonContainer();

	// copy mesh instance polygons into one big polygon soup
	mStat.polys = AddToPolygonSoup(objects, *polys);

	// construct tree from polygon soup
	Construct(polys, new BoundedRayContainer());
}


void BspTree::PreprocessPolygons(PolygonContainer &polys)
{
	// preprocess: throw out polygons coincident to the view space box (not needed)
	PolygonContainer boxPolys;
	mBbox.ExtractPolys(boxPolys);
	vector<Plane3> boxPlanes;

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

	// extract planes of box
	// TODO: can be done more elegantly than here 
	// where we first extract polygons, then compute 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
	for (int i = 0; i < (int)polys.size(); ++ i)
	{
		while (!polys[i] && (i < (int)polys.size()))
		{
			swap(polys[i], polys.back());
			polys.pop_back();
		}
	}

	CLEAR_CONTAINER(boxPolys);
}


void BspTree::Construct(const RayContainer &sampleRays,
						AxisAlignedBox3 *forcedBoundingBox)
{
	// generate new view cells for this contruction type
	mUsePredefinedViewCells = false;

    mStat.nodes = 1;

	if (forcedBoundingBox)
	{
		mBbox = *forcedBoundingBox;
	}
	else
	{
		mBbox.Initialize(); // initialise BSP tree bounding box	
	}

	PolygonContainer *polys = new PolygonContainer();
	BoundedRayContainer *rays = new BoundedRayContainer();

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

	long startTime = GetTime();
	Debug << "**** Extracting polygons from rays ****\n";

	std::map<Face *, Polygon3 *> facePolyMap;

	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		//////////////
		//-- extract polygons intersected by the rays
		Ray *ray = *rit;
	
		// get ray-face intersection. Store polygon representing the rays together
		// with rays intersecting the face.
		if (!ray->intersections.empty())
		{
			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->intersections[0].mObject);

			Mesh *mesh;
			if (obj->Type() == Intersectable::TRANSFORMED_MESH_INSTANCE)
			{
				TransformedMeshInstance *tmobj = 
					dynamic_cast<TransformedMeshInstance *>(obj);
			
				mesh = new Mesh();
				tmobj->GetTransformedMesh(*mesh);
			}
			else // MeshInstance
			{
				mesh = obj->GetMesh();
			}

			Face *face = mesh->mFaces[ray->intersections[0].mFace];
			std::map<Face *, Polygon3 *>::iterator it = facePolyMap.find(face);

			if (it != facePolyMap.end()) 
			{
				//store rays if needed for heuristics
				if (mSplitPlaneStrategy & BLOCKED_RAYS)
					(*it).second->mPiercingRays.push_back(ray);
			} 
			else
			{	
				//store rays if needed for heuristics
				Polygon3 *poly = new Polygon3(face, mesh);
				poly->mParent = obj;
				polys->push_back(poly);

				if (mSplitPlaneStrategy & BLOCKED_RAYS)
					poly->mPiercingRays.push_back(ray);

				facePolyMap[face] = poly;
			}

			// cleanup
			if (obj->Type() == Intersectable::TRANSFORMED_MESH_INSTANCE)
				DEL_PTR(mesh);
		}
	}
	// clear helper structure
	// note: memory will not be released using clear!
	facePolyMap.clear();

	// compute bounding box
	if (!forcedBoundingBox)
	{
		mBbox.Include(*polys);
	}

	////////////
	//-- store rays
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		Ray *ray = *rit;
        ray->SetId(-1); // reset id

		float minT, maxT;
		if (mBbox.GetRaySegment(*ray, minT, maxT))
			rays->push_back(new BoundedRay(ray, minT, maxT));
	}

	// throw out bad polygons
	PreprocessPolygons(*polys);

	mStat.polys = (int)polys->size();

	Debug << "**** Finished polygon extraction ****" << endl;
	Debug << (int)polys->size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
	Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;

	Construct(polys, rays);
}


void BspTree::Construct(const ObjectContainer &objects, 
						const RayContainer &sampleRays,
						AxisAlignedBox3 *forcedBoundingBox)
{
	// generate new view cells for this construction type
	mUsePredefinedViewCells = false;

    mStat.nodes = 1;
	mBbox.Initialize(); 	// initialise BSP tree bounding box
	
	if (forcedBoundingBox)
	{
		mBbox = *forcedBoundingBox;
	}

	BoundedRayContainer *rays = new BoundedRayContainer();
	PolygonContainer *polys = new PolygonContainer();
	
	// copy mesh instance polygons into one big polygon soup
	mStat.polys = AddToPolygonSoup(objects, *polys, 0, !forcedBoundingBox);

	///////
	//-- store rays
	RayContainer::const_iterator rit, rit_end = sampleRays.end();
	
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		Ray *ray = *rit;
        ray->SetId(-1); // reset id

		float minT, maxT;
		if (mBbox.GetRaySegment(*ray, minT, maxT))
			rays->push_back(new BoundedRay(ray, minT, maxT));
	}

	PreprocessPolygons(*polys);
	Debug << "tree has " << (int)polys->size() << " polys, " << (int)sampleRays.size() << " rays" << endl;
	Construct(polys, rays);
}


void BspTree::Construct(PolygonContainer *polys, BoundedRayContainer *rays)
{
	BspTraversalStack tStack;
	mRoot = new BspLeaf();

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

	BspTraversalData tData(mRoot, 
						   polys, 
						   0, 
						   mOutOfBoundsCell, 
						   rays, 
						   ComputePvsSize(*rays), 
						   mUseAreaForPvs ? geom->GetArea() : geom->GetVolume(), 
						   geom);

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

	tStack.push(tData);

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

	mStat.Start();
	cout << "Constructing bsp tree ...\n";
	long startTime = GetTime();
	while (!tStack.empty()) 
	{
		tData = tStack.top();
	    tStack.pop();

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

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

		if (mStat.Leaves() >= nleaves)
		{
			nleaves += 500;
			
			cout << "leaves=" << mStat.Leaves() << endl;
			Debug << "needed "
				  << TimeDiff(interTime, GetTime())*1e-3 
				  << " secs to create 500 leaves" << endl;
			interTime = GetTime();
		}
	}

	cout << "finished\n";

	mStat.Stop();
}


bool BspTree::TerminationCriteriaMet(const BspTraversalData &data) const
{
	return 
		(((int)data.mPolygons->size() <= mTermMinPolys) || 
		 ((int)data.mRays->size() <= mTermMinRays) ||
		 (data.mPvs <= mTermMinPvs) ||
		 (data.mProbability <= mTermMinProbability) ||
		 (data.mDepth >= mTermMaxDepth) ||
		 (mStat.Leaves() >= mMaxViewCells) ||
		 (data.GetAvgRayContribution() > mTermMaxRayContribution));
}


BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
{
	////////
	//-- terminate traversal

	if (TerminationCriteriaMet(tData))		
	{
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
	
		BspViewCell *viewCell;

		if (!mUsePredefinedViewCells)
		{	// generate new view cell for each leaf
			viewCell = new BspViewCell();cout << "g";
		}
		else
		{	
			// add predefined view cell to leaf
			viewCell = dynamic_cast<BspViewCell *>(tData.mViewCell);

			// out of bounds cell can be handled as any other view cell,
			// responsibility for deleting it has the view cells manager.
			if (IsOutOfBounds(viewCell))
			{
				mOutOfBoundsCellPartOfTree = true;
			}
		
		}

		leaf->SetViewCell(viewCell);
		viewCell->mLeaves.push_back(leaf);

		//float probability = max(0.0f, tData.mProbability);
		float probability = tData.mProbability;

		if (mUseAreaForPvs)
			viewCell->SetArea(probability);
		else
			viewCell->SetVolume(probability);
		
		///////////
		//-- add pvs contribution of rays

		if (viewCell != mOutOfBoundsCell)
		{
			int conSamp = 0, sampCon = 0;
			AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
			
			mStat.contributingSamples += conSamp;
			mStat.sampleContributions += sampCon;
		}

		if (1) EvaluateLeafStats(tData);
		
		////////
		//-- clean up

		// discard polygons
		CLEAR_CONTAINER(*tData.mPolygons);
		// discard rays
		CLEAR_CONTAINER(*tData.mRays);

		delete tData.mPolygons;
		delete tData.mRays;
		delete tData.mGeometry;

		return leaf;
	}

	///////////
	//-- continue subdivision

	PolygonContainer coincident;
	
	BspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mOutOfBoundsCell, 
								new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
	BspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mOutOfBoundsCell, 
							   new BoundedRayContainer(), 0, 0, new BspNodeGeometry());


	// create new interior node and two leaf nodes
	BspInterior *interior = 
		SubdivideNode(tData, tFrontData, tBackData, coincident);

	if (1)
	{
		int pvsData = tData.mPvs;

		float cData = (float)pvsData * tData.mProbability;
		float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
		float cBack = (float)tBackData.mPvs * tBackData.mProbability;

		float costDecr = (cFront + cBack - cData) / mBbox.GetVolume();
		
		mTotalCost += costDecr;
		mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - pvsData;

		mSubdivisionStats 
			<< "#ViewCells\n" << mStat.Leaves() << endl
			<< "#RenderCostDecrease\n" << -costDecr << endl 
			<< "#TotalRenderCost\n" << mTotalCost << endl
			<< "#AvgRenderCost\n" << mTotalPvsSize / mStat.Leaves() << endl;
	}

	// extract view cells from coincident polygons 
	// with respect to the orientation of their normal
    // note: if front or back polygons are empty, 
	// we get the valid in - out classification for the view cell

	if (mUsePredefinedViewCells)
	{
   		ExtractViewCells(tFrontData,
						 tBackData,
						 coincident,
						 interior->mPlane);			 
	}

	// don't need coincident polygons anymory
	CLEAR_CONTAINER(coincident);

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

	////////
	//-- cleanup

	DEL_PTR(tData.mNode);
	DEL_PTR(tData.mPolygons);
	DEL_PTR(tData.mRays);
	DEL_PTR(tData.mGeometry);		
	
	return interior;
}


void BspTree::ExtractViewCells(BspTraversalData &frontData,
							   BspTraversalData &backData, 
							   const PolygonContainer &coincident,
							   const Plane3 &splitPlane) const
{
	// if not empty, tree is further subdivided => don't have to find view cell
	bool foundFront = !frontData.mPolygons->empty();
	bool foundBack = !frontData.mPolygons->empty();

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

	//////////
	//-- find first view cells in front and back leafs

	for (; !(foundFront && foundBack) && (it != it_end); ++ it)
	{
		if (DotProd((*it)->GetNormal(), splitPlane.mNormal) > 0)
		{
			backData.mViewCell = dynamic_cast<ViewCellLeaf *>((*it)->mParent);
			foundBack = true;
		}
		else
		{
			frontData.mViewCell = dynamic_cast<ViewCellLeaf *>((*it)->mParent);
			foundFront = true;
		}
	}
}


BspInterior *BspTree::SubdivideNode(BspTraversalData &tData,
									BspTraversalData &frontData,
									BspTraversalData &backData,
									PolygonContainer &coincident)
{
	mStat.nodes += 2;
	
	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);	
	
	// select subdivision plane
	BspInterior *interior = new BspInterior(SelectPlane(leaf, tData)); 
	
#ifdef GTP_DEBUG
	Debug << interior << endl;
#endif
	
	// subdivide rays into front and back rays
	SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);

	// subdivide polygons with plane
	mStat.polySplits += SplitPolygons(interior->GetPlane(),
									  *tData.mPolygons, 
									  *frontData.mPolygons, 
									  *backData.mPolygons, 
									  coincident);
	

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

	// split geometry and compute area
	if (1)
	{
		tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
									   *backData.mGeometry, 
									   interior->mPlane,
									   mBbox,
									   //0.000000000001);
									   mEpsilon);
	
		
		if (mUseAreaForPvs)
		{
			frontData.mProbability = frontData.mGeometry->GetArea();
			backData.mProbability = backData.mGeometry->GetArea();
		}
		else
		{
			frontData.mProbability = frontData.mGeometry->GetVolume();
			backData.mProbability = tData.mProbability - frontData.mProbability;
		}
	}

	//-- create front and back leaf

	BspInterior *parent = leaf->GetParent();

	// replace a link from node's parent
	if (!leaf->IsRoot())
	{
		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 BspTree::SortSubdivisionCandidates(const PolygonContainer &polys, 
								  const int axis, 
								  vector<SortableEntry> &splitCandidates) const
{
	splitCandidates.clear();

	int requestedSize = 2 * (int)polys.size();
	// creates a sorted split candidates array  
	splitCandidates.reserve(requestedSize);

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

	AxisAlignedBox3 box;

	// insert all queries 
	for(it = polys.begin(); it != it_end; ++ it)
	{
		box.Initialize();
		box.Include(*(*it));
		
		splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
		splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
	}

	stable_sort(splitCandidates.begin(), splitCandidates.end());
}


float BspTree::BestCostRatio(const PolygonContainer &polys,
							 const AxisAlignedBox3 &box,
							 const int axis,
							 float &position,
							 int &objectsBack,
							 int &objectsFront) const
{
	vector<SortableEntry> splitCandidates;

	SortSubdivisionCandidates(polys, axis, splitCandidates);
	
	// go through the lists, count the number of objects left and right
	// and evaluate the following cost funcion:
	// C = ct_div_ci  + (ol + or)/queries
	
	int objectsLeft = 0, objectsRight = (int)polys.size();
	
	float minBox = box.Min(axis);
	float maxBox = box.Max(axis);
	float boxArea = box.SurfaceArea();
  
	float minBand = minBox + mSplitBorder * (maxBox - minBox);
	float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
	
	float minSum = 1e20f;
	vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();

	for(ci = splitCandidates.begin(); ci != ci_end; ++ ci) 
	{
		switch ((*ci).type) 
		{
			case SortableEntry::POLY_MIN:
				++ objectsLeft;
				break;
			case SortableEntry::POLY_MAX:
			    -- objectsRight;
				break;
			default:
				break;
		}
		
		if ((*ci).value > minBand && (*ci).value < maxBand) 
		{
			AxisAlignedBox3 lbox = box;
			AxisAlignedBox3 rbox = box;
			lbox.SetMax(axis, (*ci).value);
			rbox.SetMin(axis, (*ci).value);

			const float sum = objectsLeft * lbox.SurfaceArea() + 
					 		  objectsRight * rbox.SurfaceArea();
      
			if (sum < minSum) 
			{
				minSum = sum;
				position = (*ci).value;

				objectsBack = objectsLeft;
				objectsFront = objectsRight;
			}
		}
	}
  
	const float oldCost = (float)polys.size();
	const float newCost = mAxisAlignedCtDivCi + minSum / boxArea;
	const float ratio = newCost / oldCost;

#ifdef GTP_DEBUG
  Debug << "====================" << endl;
  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
#endif
  return ratio;
}

bool BspTree::SelectAxisAlignedPlane(Plane3 &plane, 
									 const PolygonContainer &polys) const
{
	AxisAlignedBox3 box;
	box.Initialize();
	
	// create bounding box of region
	box.Include(polys);
        
	int objectsBack = 0, objectsFront = 0;
	int axis = 0;
	float costRatio = MAX_FLOAT;
	Vector3 position;

	//-- area subdivision
	for (int i = 0; i < 3; ++ i) 
	{
		float p = 0;
		float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
		
		if (r < costRatio)
		{
			costRatio = r;
			axis = i;
			position = p;
		}
	}
	
	if (costRatio >= mMaxCostRatio)
		return false;

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

	return true;
}


Plane3 BspTree::SelectPlane(BspLeaf *leaf, BspTraversalData &data)
{
	if ((!mMaxPolyCandidates || data.mPolygons->empty()) && 
		(!mMaxRayCandidates || data.mRays->empty()))
	{
		Debug << "Warning: No autopartition polygon candidate available\n";
	
		// return axis aligned split
		AxisAlignedBox3 box;
		box.Initialize();
	
		// create bounding box of region
		box.Include(*data.mPolygons);

		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;
		return Plane3(norm, position);
	}
	
	if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
		((int)data.mPolygons->size() > mTermMinPolysForAxisAligned) &&
		((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
		((mTermMinObjectsForAxisAligned < 0) || 
		 (Polygon3::ParentObjectsSize(*data.mPolygons) > mTermMinObjectsForAxisAligned)))
	{
		Plane3 plane;
		if (SelectAxisAlignedPlane(plane, *data.mPolygons))
			return plane;
	}

	// simplest strategy: just take next polygon
	if (mSplitPlaneStrategy & RANDOM_POLYGON)
	{
        if (!data.mPolygons->empty())
		{
			Polygon3 *nextPoly = 
				(*data.mPolygons)[(int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1))];
			return nextPoly->GetSupportingPlane();
		}
		else
		{
			const int candidateIdx = (int)RandomValue(0, (Real)((int)data.mRays->size() - 1));
			BoundedRay *bRay = (*data.mRays)[candidateIdx];

			Ray *ray = bRay->mRay;
						
			const Vector3 minPt = ray->Extrap(bRay->mMinT);
			const Vector3 maxPt = ray->Extrap(bRay->mMaxT);

			const Vector3 pt = (maxPt + minPt) * 0.5;

			const Vector3 normal = ray->GetDir();
			
			return Plane3(normal, pt);
		}

		return Plane3();
	}

	// use heuristics to find appropriate plane
	return SelectPlaneHeuristics(leaf, data); 
}


Plane3 BspTree::ChooseCandidatePlane(const BoundedRayContainer &rays) const
{	
	const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
	BoundedRay *bRay = rays[candidateIdx];
	Ray *ray = bRay->mRay;

	const Vector3 minPt = ray->Extrap(bRay->mMinT);
	const Vector3 maxPt = ray->Extrap(bRay->mMaxT);

	const Vector3 pt = (maxPt + minPt) * 0.5;

	const Vector3 normal = ray->GetDir();
			
	return Plane3(normal, pt);
}

Plane3 BspTree::ChooseCandidatePlane2(const BoundedRayContainer &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);

		BoundedRay *bRay = rays[idx[j]];
		pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : 
							bRay->mRay->Extrap(bRay->mMaxT);
	}	
			
	return Plane3(pt[0], pt[1], pt[2]);
}


Plane3 BspTree::ChooseCandidatePlane3(const BoundedRayContainer &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 BoundedRay *ray1 = rays[idx1];
	const BoundedRay *ray2 = rays[idx2];

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

	const Vector3 orig1 = ray1->mRay->Extrap(ray1->mMinT);
	const Vector3 orig2 = ray2->mRay->Extrap(ray2->mMinT);

	// vector from line 1 to line 2
	const Vector3 vd = orig1 - orig2;
	
	// 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 = orig1 + norm * dist * 0.5;

	return Plane3(norm, planePt);
}


Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf, BspTraversalData &data)
{
	float lowestCost = MAX_FLOAT;
	Plane3 bestPlane;
	// intermediate plane
	Plane3 plane;

	const int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
	int maxIdx = (int)data.mPolygons->size();
	
	for (int i = 0; i < limit; ++ i)
	{
		// 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
		const float candidateCost = 
			SplitPlaneCost(poly->GetSupportingPlane(), data);

		if (candidateCost < lowestCost)
		{
			bestPlane = poly->GetSupportingPlane();
			lowestCost = candidateCost;
		}
	}
	
	//-- choose candidate planes extracted from rays
	for (int i = 0; i < mMaxRayCandidates; ++ i)
	{
		plane = ChooseCandidatePlane3(*data.mRays);
		const float candidateCost = SplitPlaneCost(plane, data);

		if (candidateCost < lowestCost)
		{
			bestPlane = plane;	
			lowestCost = candidateCost;
		}
	}

#ifdef GTP_DEBUG
	Debug << "plane lowest cost: " << lowestCost << endl;
#endif
	
	return bestPlane;
}


float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
							  const PolygonContainer &polys) const
{
	float val = 0;

	float sumBalancedPolys = 0;
	float sumSplits = 0;
	float sumPolyArea = 0;
	float sumBalancedViewCells = 0;
	float sumBlockedRays = 0;
	float totalBlockedRays = 0;
	//float totalArea = 0;
	int totalViewCells = 0;

	// need three unique ids for each type of view cell
	// for balanced view cells criterium
	ViewCell::NewMail();
	const int backId = ViewCell::sMailId;
	ViewCell::NewMail();
	const int frontId = ViewCell::sMailId;
	ViewCell::NewMail();
	const int frontAndBackId = ViewCell::sMailId;

	bool useRand;
	int limit;

	// choose test polyongs randomly if over threshold
	if ((int)polys.size() > mMaxTests)
	{
		useRand = true;
		limit = mMaxTests;
	}
	else
	{
		useRand = false;
		limit = (int)polys.size();
	}

	for (int i = 0; i < limit; ++ i)
	{
		const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;

		Polygon3 *poly = polys[testIdx];

        const int classification = 
			poly->ClassifyPlane(candidatePlane, mEpsilon);

		if (mSplitPlaneStrategy & BALANCED_POLYS)
			sumBalancedPolys += sBalancedPolysTable[classification];
		
		if (mSplitPlaneStrategy & LEAST_SPLITS)
			sumSplits += sLeastPolySplitsTable[classification];

		if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
		{
			if (classification == Polygon3::COINCIDENT)
				sumPolyArea += poly->GetArea();
			//totalArea += area;
		}
		
		if (mSplitPlaneStrategy & BLOCKED_RAYS)
		{
			const float blockedRays = (float)poly->mPiercingRays.size();
		
			if (classification == Polygon3::COINCIDENT)
				sumBlockedRays += blockedRays;
			
			totalBlockedRays += blockedRays;
		}

		// assign view cells to back or front according to classificaion
		if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
		{
			MeshInstance *viewCell = poly->mParent;
		
			// assure that we only count a view cell 
			// once for the front and once for the back side of the plane
			if (classification == Polygon3::FRONT_SIDE)
			{
				if ((viewCell->mMailbox != frontId) && 
					(viewCell->mMailbox != frontAndBackId))
				{
					sumBalancedViewCells += 1.0;

					if (viewCell->mMailbox != backId)
						viewCell->mMailbox = frontId;
					else
						viewCell->mMailbox = frontAndBackId;
					
					++ totalViewCells;
				}
			}
			else if (classification == Polygon3::BACK_SIDE)
			{
				if ((viewCell->mMailbox != backId) &&
				    (viewCell->mMailbox != frontAndBackId))
				{
					sumBalancedViewCells -= 1.0;

					if (viewCell->mMailbox != frontId)
						viewCell->mMailbox = backId;
					else
						viewCell->mMailbox = frontAndBackId; 

					++ totalViewCells;
				}
			}
		}
	}

	const float polysSize = (float)polys.size() + Limits::Small;

	// all values should be approx. between 0 and 1 so they can be combined
	// and scaled with the factors according to their importance
	if (mSplitPlaneStrategy & BALANCED_POLYS)
		val += mBalancedPolysFactor * fabs(sumBalancedPolys) / polysSize;
	
	if (mSplitPlaneStrategy & LEAST_SPLITS)  
		val += mLeastSplitsFactor * sumSplits / polysSize;

	if (mSplitPlaneStrategy & LARGEST_POLY_AREA) 
		// HACK: polys.size should be total area so scaling is between 0 and 1
		val += mLargestPolyAreaFactor * (float)polys.size() / sumPolyArea;

	if (mSplitPlaneStrategy & BLOCKED_RAYS)
		if (totalBlockedRays != 0)
			val += mBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;

	if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
		val += mBalancedViewCellsFactor * fabs(sumBalancedViewCells) / 
			((float)totalViewCells + Limits::Small);
	
	return val;
}


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


float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
							  const BoundedRayContainer &rays,
							  const int pvs,
							  const float probability,
							  const BspNodeGeometry &cell) const
{
	float val = 0;

	float sumBalancedRays = 0;
	float sumRaySplits = 0;

	int frontPvs = 0;
	int backPvs = 0;

	// probability that view point lies in child
	float pOverall = 0;
	float pFront = 0;
	float pBack = 0;

	const bool pvsUseLen = false;

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

		// construct child geometry with regard to the candidate split plane
		BspNodeGeometry geomFront;
		BspNodeGeometry geomBack;

		const bool splitSuccessFull = 
			cell.SplitGeometry(geomFront, 
							   geomBack, 
							   candidatePlane,
							   mBbox,
							   mEpsilon);

		if (mUseAreaForPvs)
		{
			pFront = geomFront.GetArea();
			pBack = geomBack.GetArea();
		}
		else
		{
			pFront = geomFront.GetVolume();
			pBack = pOverall - pFront;
		}
		
		
		// give penalty to unbalanced split
		if (1 && 
			(!splitSuccessFull || (pFront <= 0) || (pBack <= 0) || 
			!geomFront.Valid() || !geomBack.Valid()))
		{
			//Debug << "error f: " << pFront << " b: " << pBack << endl;
			return 99999.9f;
		}

		pOverall = probability;
	}
			
	bool useRand;
	int limit;

	// choose test polyongs randomly if over threshold
	if ((int)rays.size() > mMaxTests)
	{
		useRand = true;
		limit = mMaxTests;
	}
	else
	{
		useRand = false;
		limit = (int)rays.size();
	}

	for (int i = 0; i < limit; ++ i)
	{
		const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;
	
		BoundedRay *bRay = rays[testIdx];

		Ray *ray = bRay->mRay;
		const float minT = bRay->mMinT;
		const float maxT = bRay->mMaxT;

		Vector3 entP, extP;

		const int cf = 
			ray->ClassifyPlane(candidatePlane, minT, maxT, entP, extP);

		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
		{
			sumBalancedRays += sBalancedRaysTable[cf];
		}
		
		if (mSplitPlaneStrategy & BALANCED_RAYS)
		{
			sumRaySplits += sLeastRaySplitsTable[cf];
		}

		if (mSplitPlaneStrategy & PVS)
		{
			// in case the ray intersects an object
			// assure that we only count the object 
			// once for the front and once for the back side of the plane
			
			// add the termination object
			if (!ray->intersections.empty())
				AddObjToPvs(ray->intersections[0].mObject, cf, frontPvs, backPvs);
			
			// add the source object
			AddObjToPvs(ray->sourceObject.mObject, cf, frontPvs, backPvs);
		}
	}

	const float raysSize = (float)rays.size() + Limits::Small;

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

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

	const float denom = pOverall * (float)pvs + Limits::Small;

	if (mSplitPlaneStrategy & PVS)
	{
		val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
	}

	
#ifdef GTP_DEBUG
	Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
		  << " frontpvs: " << frontPvs << " pFront: " << pFront 
		  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
#endif
	
	return val;
}

void BspTree::AddObjToPvs(Intersectable *obj,
						  const int cf,
						  int &frontPvs,
						  int &backPvs) const
{
	if (!obj)
		return;
	// TODO: does this really belong to no pvs?
	//if (cf == Ray::COINCIDENT) return;

	// object belongs to both PVS
	const bool bothSides = (cf == Ray::FRONT_BACK) || 
						   (cf == Ray::BACK_FRONT) ||
						   (cf == Ray::COINCIDENT);

	if ((cf == Ray::FRONT) || bothSides)
	{
		if ((obj->mMailbox != sFrontId) && 
			(obj->mMailbox != sFrontAndBackId))
		{
			++ frontPvs;

			if (obj->mMailbox == sBackId)
				obj->mMailbox = sFrontAndBackId;	
			else
				obj->mMailbox = sFrontId;								
		}
	}
	
	if ((cf == Ray::BACK) || bothSides)
	{
		if ((obj->mMailbox != sBackId) &&
			(obj->mMailbox != sFrontAndBackId))
		{
			++ backPvs;

			if (obj->mMailbox == sFrontId)
				obj->mMailbox = sFrontAndBackId; 
			else
				obj->mMailbox = sBackId;				
		}
	}
}


float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
							  BspTraversalData &data) const
{
	float val = 0;

	if (mSplitPlaneStrategy & VERTICAL_AXIS)
	{
		Vector3 tinyAxis(0,0,0); tinyAxis[mBbox.Size().TinyAxis()] = 1.0f;
		// we put a penalty on the dot product between the "tiny" vertical axis
		// and the split plane axis
		val += mVerticalSplitsFactor * 
			   fabs(DotProd(candidatePlane.mNormal, tinyAxis));
	}

	// the following criteria loop over all polygons to find the cost value
	if ((mSplitPlaneStrategy & BALANCED_POLYS)      ||
		(mSplitPlaneStrategy & LEAST_SPLITS)        ||
		(mSplitPlaneStrategy & LARGEST_POLY_AREA)   ||
		(mSplitPlaneStrategy & BALANCED_VIEW_CELLS) ||
		(mSplitPlaneStrategy & BLOCKED_RAYS))
	{
		val += SplitPlaneCost(candidatePlane, *data.mPolygons);
	}

	// the following criteria loop over all rays to find the cost value
	if ((mSplitPlaneStrategy & BALANCED_RAYS)      ||
		(mSplitPlaneStrategy & LEAST_RAY_SPLITS)   ||
		(mSplitPlaneStrategy & PVS))
	{
		val += SplitPlaneCost(candidatePlane, *data.mRays, data.mPvs, 
							  data.mProbability, *data.mGeometry);
	}

	// return linear combination of the sums
	return val;
}

void BspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
{
	stack<BspNode *> nodeStack;
	nodeStack.push(mRoot);
 
	while (!nodeStack.empty()) 
	{
		BspNode *node = nodeStack.top();
    
		nodeStack.pop();
    
		if (node->IsLeaf()) 
		{
			BspLeaf *leaf = (BspLeaf *)node;		
			leaves.push_back(leaf);
		} 
		else 
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

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


AxisAlignedBox3 BspTree::GetBoundingBox() const
{
	return mBbox;
}


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

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

	if (data.mDepth < mStat.minDepth)
		mStat.minDepth = data.mDepth;

	// accumulate depth to compute average depth
	mStat.accumDepth += data.mDepth;
	// accumulate rays to compute rays /  leaf
	mStat.accumRays += (int)data.mRays->size();

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

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

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

	if (data.GetAvgRayContribution() > mTermMaxRayContribution)
		++ mStat.maxRayContribNodes;
	
	if (data.mProbability <= mTermMinProbability) 
		++ mStat.minProbabilityNodes;

#ifdef GTP_DEBUG
	Debug << "BSP stats: "
		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
		  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
		  << "Probability: " << data.mProbability << " (min: " << mTermMinProbability << "), "
		  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << mTermMinPolys << "), "
		  << "#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 BspTree::_CastRay(Ray &ray)
{
	int hits = 0;
  
	stack<BspRayTraversalData> tStack;
  
	float maxt, mint;

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

	ViewCell::NewMail();

	Vector3 entp = ray.Extrap(mint);
	Vector3 extp = ray.Extrap(maxt);
  
	BspNode *node = mRoot;
	BspNode *farChild = NULL;
	
	while (1)
	{
		if (!node->IsLeaf()) 
		{
			BspInterior *in = dynamic_cast<BspInterior *>(node);
			
			Plane3 splitPlane = in->GetPlane();
			const int entSide = splitPlane.Side(entp);
			const int extSide = splitPlane.Side(extp);

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

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

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

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

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

			// push data for far child
			tStack.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->mViewCell->Mailed())
			{
				//ray.bspIntersections.push_back(Ray::BspIntersection(maxt, leaf));
				leaf->mViewCell->Mail();
				++ hits;
			}
			
   			//-- fetch the next far child from the stack
			if (tStack.empty())
				break;
      
			entp = extp;
			mint = maxt; // NOTE: need this?

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

			BspRayTraversalData &s = tStack.top();

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

			tStack.pop();
		}
	}

	return hits;
}


int BspTree::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;


	const float thresh = 1 ? 1e-6f : 0.0f;

	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, maxt));
		
			// find intersection of ray segment with plane
			float t;
			extp = splitPlane.FindIntersection(origin, extp, &t);
			maxt *= t;  
		} 
		else 
		{
			// reached leaf => intersection with view cell
			BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
		
			if (!leaf->mViewCell->Mailed()) 
			{
				viewcells.push_back(leaf->mViewCell);
				leaf->mViewCell->Mail();
				++ hits;
			}
		
			//-- fetch the next far child from the stack
			if (tStack.empty())
				break;
	    
			entp = extp;
			mint = maxt; // NOTE: need this?
		
			BspRayTraversalData &s = tStack.top();
		
			node = s.mNode;
			extp = s.mExitPoint;
			maxt = s.mMaxT;
		
			tStack.pop();
		}
	}

	return hits;
}


void BspTree::CollectViewCells(ViewCellContainer &viewCells) const
{
	stack<BspNode *> nodeStack;
	nodeStack.push(mRoot);

	ViewCell::NewMail();

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

		if (node->IsLeaf()) 
		{
			ViewCellLeaf *viewCell = dynamic_cast<BspLeaf *>(node)->mViewCell;

			if (!viewCell->Mailed()) 
			{
				viewCell->Mail();
				viewCells.push_back(viewCell);
			}
		}
		else 
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

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


BspTreeStatistics &BspTree::GetStat()
{
	return mStat;
}


float BspTree::AccumulatedRayLength(BoundedRayContainer &rays) const
{
	float len = 0;

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

	for (it = rays.begin(); it != it_end; ++ it)
	{
		len += SqrDistance((*it)->mRay->Extrap((*it)->mMinT), 
						   (*it)->mRay->Extrap((*it)->mMaxT));
	}

	return len;
}


int BspTree::SplitRays(const Plane3 &plane,
					   BoundedRayContainer &rays, 
					   BoundedRayContainer &frontRays, 
					   BoundedRayContainer &backRays)
{
	int splits = 0;
	
	while (!rays.empty())
	{
		BoundedRay *bRay = rays.back();
		Ray *ray = bRay->mRay;
		float minT = bRay->mMinT;
		float maxT = bRay->mMaxT;

		rays.pop_back();
	
		Vector3 entP, extP;

		const int cf = 
			ray->ClassifyPlane(plane, minT, maxT, entP, extP);
		
		// set id to ray classification
		ray->SetId(cf);

		switch (cf)
		{
		case Ray::COINCIDENT: // TODO: should really discard ray?
			frontRays.push_back(bRay);
			break;
		case Ray::BACK:
			backRays.push_back(bRay);
			break;
		case Ray::FRONT:
			frontRays.push_back(bRay);
			break;
		case Ray::FRONT_BACK:
			{
				// find intersection of ray segment with plane
				const float t = plane.FindT(ray->GetLoc(), extP);
				
				const float newT = t * maxT;

				frontRays.push_back(new BoundedRay(ray, minT, newT));
				backRays.push_back(new BoundedRay(ray, newT, maxT));

				DEL_PTR(bRay);
			}
			break;
		case Ray::BACK_FRONT:
			{
				// find intersection of ray segment with plane
				const float t = plane.FindT(ray->GetLoc(), extP);
				const float newT = t * bRay->mMaxT;

				backRays.push_back(new BoundedRay(ray, minT, newT));
				frontRays.push_back(new BoundedRay(ray, newT, maxT));

				DEL_PTR(bRay);

				++ splits;
			}
			break;
		default:
			Debug << "Should not come here 4" << endl;
			break;
		}
	}

	return splits;
}

void BspTree::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 BspTree::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)
	{
		// per definition out of bounds cell has zero volume
		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);
		}
	}
}


void BspTree::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(mBbox.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]);
		}
	}
}


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


typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;


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


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

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

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

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

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

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

						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(),
												   mEpsilon);
			
			BspNode *front = interior->GetFront();
			BspNode *back = interior->GetBack();
            
			BspNodeGeometry *fGeom = new BspNodeGeometry();
			BspNodeGeometry *bGeom = new BspNodeGeometry();

			geom->SplitGeometry(*fGeom,
								*bGeom,
								interior->GetPlane(),
								mBbox,
								//0.0000001f);
								mEpsilon);
		
			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 *BspTree::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;
}


void BspTree::AddToPvs(BspLeaf *leaf,
					   const BoundedRayContainer &rays, 
					   int &sampleContributions,
					   int &contributingSamples)
{
	sampleContributions = 0;
	contributingSamples = 0;

    BoundedRayContainer::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)
	{
		int contribution = 0;
		Ray *ray = (*it)->mRay;
		float relContribution;
		if (!ray->intersections.empty())
		  contribution += vc->GetPvs().AddSample(ray->intersections[0].mObject,
												 1.0f);
												 //relContribution);
		
		if (ray->sourceObject.mObject)
			contribution += vc->GetPvs().AddSample(ray->sourceObject.mObject,
												   1.0f);
												   //relContribution);
		
		if (contribution)
		{
			sampleContributions += contribution;
			++ contributingSamples;
		}

		//if (ray->mFlags & Ray::STORE_BSP_INTERSECTIONS)
		//	ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
	}
}


int BspTree::ComputePvsSize(const BoundedRayContainer &rays) const
{
	int pvsSize = 0;

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

	Intersectable::NewMail();

	for (rit = rays.begin(); rit != rays.end(); ++ rit)
	{
		Ray *ray = (*rit)->mRay;
		
		if (!ray->intersections.empty())
		{
			if (!ray->intersections[0].mObject->Mailed())
			{
				ray->intersections[0].mObject->Mail();
				++ pvsSize;
			}
		}
		if (ray->sourceObject.mObject)
		{
			if (!ray->sourceObject.mObject->Mailed())
			{
				ray->sourceObject.mObject->Mail();
				++ pvsSize;
			}
		}
	}

	return pvsSize;
}


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


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

	int numCandidates = 0;

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

		vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();

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

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

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

	return (int)leaves.size();
}


int BspTree::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];

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

#if 0
	//////////
	//-- collect the leaves which haven't been found by ray casting
	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);
#endif

	return numLeaves;
}



/***************************************************************/
/*              BspNodeGeometry Implementation                 */
/***************************************************************/


BspNodeGeometry::BspNodeGeometry(const BspNodeGeometry &rhs)
{
	mPolys.reserve(rhs.mPolys.size());
	mPlanes.reserve(rhs.mPolys.size());

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

	int i = 0;

	for (it = rhs.mPolys.begin(); it != it_end; ++ it, ++i)
	{
		Polygon3 *poly = *it;
		Add(new Polygon3(*poly), rhs.mPlanes[i]);
	}
}


BspNodeGeometry& BspNodeGeometry::operator=(const BspNodeGeometry& g) 
{
    if (this == &g) 
		return *this;
 
	CLEAR_CONTAINER(mPolys);

	mPolys.reserve(g.mPolys.size());
	mPlanes.reserve(g.mPolys.size());

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

	int i = 0;

	for (it = g.mPolys.begin(); it != it_end; ++ it, ++ i)
	{
		Polygon3 *poly = *it;
		Add(new Polygon3(*poly), g.mPlanes[i]);
	}

	return *this;
}


BspNodeGeometry::BspNodeGeometry(const PolygonContainer &polys)
{
	mPolys = polys;
	mPlanes.reserve(polys.size());
	
	PolygonContainer::const_iterator it, it_end = polys.end();

	for (it = polys.begin(); it != it_end; ++ it)
	{
		Polygon3 *poly = *it;
		mPlanes.push_back(poly->GetSupportingPlane());
	}
}


BspNodeGeometry::~BspNodeGeometry()
{
	CLEAR_CONTAINER(mPolys);
}


int BspNodeGeometry::Size() const
{
	return (int)mPolys.size();
}


const PolygonContainer &BspNodeGeometry::GetPolys()
{
	return mPolys;
}


void BspNodeGeometry::Add(Polygon3 *p, const Plane3 &plane)
{
    mPolys.push_back(p);
	mPlanes.push_back(plane);
}


float BspNodeGeometry::GetArea() const 
{ 
	return Polygon3::GetArea(mPolys);
}


float BspNodeGeometry::GetVolume() const
{
	//-- compute volume using tetrahedralization of the geometry
	//   and adding the volume of the single tetrahedrons
	float volume = 0;
	const float f = 1.0f / 6.0f;

	PolygonContainer::const_iterator pit, pit_end = mPolys.end();

	// note: can take arbitrary point, e.g., the origin. However,
	// we rather take the center of mass to prevents precision errors
	const Vector3 center = CenterOfMass();

	for (pit = mPolys.begin(); pit != pit_end; ++ pit)
	{
		Polygon3 *poly = *pit;
		const Vector3 v0 = poly->mVertices[0] - center;

		for (int i = 1; i < (int)poly->mVertices.size() - 1; ++ i)
		{
			const Vector3 v1 = poly->mVertices[i] - center;
			const Vector3 v2 = poly->mVertices[i + 1] - center;

			// more robust version using abs and the center of mass
			volume += fabs (f * (DotProd(v0, CrossProd(v1, v2))));
		}
	}

	return volume;
}


void BspNodeGeometry::GetBoundingBox(AxisAlignedBox3 &box) const
{
	box.Initialize();
	box.Include(mPolys);
}


int BspNodeGeometry::Side(const Plane3 &plane, const float eps) const
{
	PolygonContainer::const_iterator it, it_end = mPolys.end();

	bool onFrontSide = false;
	bool onBackSide = false;

	for (it = mPolys.begin(); it != it_end; ++ it)
	{
        const int side = (*it)->Side(plane, eps);

		if (side == -1)
			onBackSide = true;
		else if (side == 1)
			onFrontSide = true;
		
		if ((side == 0) || (onBackSide && onFrontSide))
			return 0;
	}

	if (onFrontSide)
		return 1;

	return -1;
}


int BspNodeGeometry::ComputeIntersection(const AxisAlignedBox3 &box) const
{
	// 1: box does not intersect geometry
	// 0: box intersects geometry
	// -1: box contains geometry

	AxisAlignedBox3 boundingBox;
	GetBoundingBox(boundingBox);

	// no intersections with bounding box
	if (!Overlap(box, boundingBox))
	{
		return 1;
	}

	// box cotains bounding box of geometry >=> contains geometry
	if (box.Includes(boundingBox))
		return -1;

	int result = 0;

	// test geometry planes for intersections
	vector<Plane3>::const_iterator it, it_end = mPlanes.end();

	for (it = mPlanes.begin(); it != it_end; ++ it)
	{		
		const int side = box.Side(*it);
		
		// box does not intersects geometry
		if (side == 1)
		{
			return 1;
		}
		//if (side == 0) result = 0;
  }

  return result;
}


Vector3 BspNodeGeometry::CenterOfMass() const
{
	int n = 0;

	Vector3 center(0,0,0);

	PolygonContainer::const_iterator pit, pit_end = mPolys.end();

	for (pit = mPolys.begin(); pit != pit_end; ++ pit)
	{
		Polygon3 *poly = *pit;
		
		VertexContainer::const_iterator vit, vit_end = poly->mVertices.end();

		for(vit = poly->mVertices.begin(); vit != vit_end; ++ vit)
		{
			center += *vit;
			++ n;
		}
	}

	//Debug << "center: " << center << " new " << center / (float)n << endl;

	return center / (float)n;
}


bool BspNodeGeometry::Valid() const
{
	// geometry is degenerated
	if (mPolys.size() < 4)
		return false;
	
	return true;
}

void IncludeNodeGeomInMesh(const BspNodeGeometry &geom, Mesh &mesh)
{
	// add single polygons to mesh
	PolygonContainer::const_iterator it, it_end = geom.mPolys.end();
	
	for (it = geom.mPolys.begin(); it != geom.mPolys.end(); ++ it)
	{
		Polygon3 *poly = (*it);
		IncludePolyInMesh(*poly, mesh);
	}
}


bool BspNodeGeometry::SplitGeometry(BspNodeGeometry &front,
									BspNodeGeometry &back,
									const Plane3 &splitPlane,
									const AxisAlignedBox3 &box,
									const float epsilon) const
{	
	//-- trivial cases
	if (0)
	{
		const int cf = Side(splitPlane, epsilon);

		if (cf == -1)
		{
			back = *this;
			return false;
		}
		else if (cf == 1)
		{
			front = *this;
			return false;
		}
	}

	// get cross section of new polygon
	Polygon3 *planePoly = box.CrossSection(splitPlane);
	
	// split polygon with all other polygons
	planePoly = SplitPolygon(planePoly, epsilon);

	bool splitsGeom = (planePoly != NULL);

	//-- new polygon splits all other polygons
	for (int i = 0; i < (int)mPolys.size()/* && planePoly*/; ++ i)
	{
		/// don't use epsilon here to get exact split planes
		const int cf = 
			mPolys[i]->ClassifyPlane(splitPlane, epsilon);
			
		switch (cf)
		{
			case Polygon3::SPLIT:
				{
					Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);

					Polygon3 *frontPoly = new Polygon3();
					Polygon3 *backPoly = new Polygon3();
				
					poly->Split(splitPlane, 
								*frontPoly, 
								*backPoly,
								epsilon);

					DEL_PTR(poly);

					if (frontPoly->Valid(epsilon))
					{
						front.Add(frontPoly, mPlanes[i]);
					}
					else
					{
						//Debug << "no f! " << endl;
						DEL_PTR(frontPoly);
					}

					if (backPoly->Valid(epsilon))
					{
						back.Add(backPoly, mPlanes[i]);
					}
					else
					{
						//Debug << "no b! " << endl;
						DEL_PTR(backPoly);
					}
				}
				
				break;
			case Polygon3::BACK_SIDE:
				back.Add(new Polygon3(mPolys[i]->mVertices), mPlanes[i]);
				break;
			case Polygon3::FRONT_SIDE:
				front.Add(new Polygon3(mPolys[i]->mVertices), mPlanes[i]);
				break;
			// only put into back container (split should have no effect ...)
			case Polygon3::COINCIDENT:
				//Debug << "error should not come here" << endl;
				splitsGeom = false;

				if (DotProd(mPlanes[i].mNormal, splitPlane.mNormal > 0))
					back.Add(new Polygon3(mPolys[i]->mVertices), mPlanes[i]);
				else
                    front.Add(new Polygon3(mPolys[i]->mVertices), mPlanes[i]);
				break;
			default:
				break;
		}
	}

	//-- finally add the new polygon to the child node geometries
	if (planePoly)
	{
        // add polygon with reverse orientation to front cell
		Plane3 reversePlane(splitPlane);
		reversePlane.ReverseOrientation();

		// add polygon with normal pointing into positive half space to back cell
		back.Add(planePoly, splitPlane);
		//back.mPolys.push_back(planePoly);

		Polygon3 *reversePoly = planePoly->CreateReversePolygon();
		front.Add(reversePoly, reversePlane);
		//Debug << "poly normal : " << reversePoly->GetSupportingPlane().mNormal << " split plane normal " << reversePlane.mNormal << endl;
	}
	

	return splitsGeom;
}


Polygon3 *BspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
										const float epsilon) const
{
	if (!planePoly->Valid(epsilon))
	{
		//Debug << "not valid!!" << endl;
		DEL_PTR(planePoly);
	}
	// polygon is split by all other planes
	for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
	{
		//Plane3 plane = mPolys[i]->GetSupportingPlane();
		Plane3 plane = mPlanes[i];

		/// don't use epsilon here to get exact split planes
		const int cf = 
			planePoly->ClassifyPlane(plane, epsilon);
			
		// split new polygon with all previous planes
		switch (cf)
		{
			case Polygon3::SPLIT:
				{
					Polygon3 *frontPoly = new Polygon3();
					Polygon3 *backPoly = new Polygon3();

					planePoly->Split(plane, 
									 *frontPoly, 
									 *backPoly,
									 epsilon);
					
					// don't need anymore
					DEL_PTR(planePoly);
					DEL_PTR(frontPoly);

					// back polygon is belonging to geometry
					if (backPoly->Valid(epsilon))
						planePoly = backPoly;
					else
						DEL_PTR(backPoly);
				}
				break;
			case Polygon3::FRONT_SIDE:
				DEL_PTR(planePoly);
                break;
			// polygon is taken as it is
			case Polygon3::BACK_SIDE:
			case Polygon3::COINCIDENT:
			default:
				break;
		}
	}

	return planePoly;
}


ViewCell *BspTree::GetViewCell(const Vector3 &point)
{
	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)->mViewCell;
			break;
		} 
		else 
		{

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

			// random decision
			if (interior->GetPlane().Side(point) < 0)
			{
				nodeStack.push(interior->GetBack());
			}
			else
			{
				nodeStack.push(interior->GetFront());
			}
		}
	}

	return viewcell;
}


bool BspTree::Export(OUT_STREAM &stream)
{
	ExportNode(mRoot, stream);
	return true;
}


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

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


}