#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 <stack>
#include <time.h>
#include <iomanip>
#include "Exporter.h"
#include "Plane3.h"

int BspTree::sTermMaxPolygons = 10;
int BspTree::sTermMinPvs = 20;
int BspTree::sTermMaxDepth = 20;
int BspTree::sMaxPolyCandidates = 10;
int BspTree::sMaxRayCandidates = 10;
int BspTree::sSplitPlaneStrategy = BALANCED_POLYS; 
int BspTree::sConstructionMethod = FROM_INPUT_VIEW_CELLS;
int BspTree::sTermMaxPolysForAxisAligned = 50;
int BspTree::sTermMaxObjectsForAxisAligned = 50;
int BspTree::sTermMaxRaysForAxisAligned = -1;
int BspTree::sTermMaxRays = -1;
int BspTree::sMinPvsDif = 10;
int BspTree::sMinPvs = 10;
int BspTree::sMaxPvs = 500;

float BspTree::sCt_div_ci = 1.0f;
float BspTree::sSplitBorder = 1.0f;
float BspTree::sMaxCostRatio = 0.9f;

//-- factors for bsp tree split plane heuristics

float BspTree::sLeastSplitsFactor = 1.0f;
float BspTree::sBalancedPolysFactor = 1.0f;
float BspTree::sBalancedViewCellsFactor = 1.0f;

// NOTE:  very important criterium for 2.5d scenes
float BspTree::sVerticalSplitsFactor = 1.0f;
float BspTree::sLargestPolyAreaFactor = 1.0f;
float BspTree::sBlockedRaysFactor = 1.0f;
float BspTree::sLeastRaySplitsFactor = 1.0f;
float BspTree::sBalancedRaysFactor = 1.0f;
float BspTree::sPvsFactor = 1.0f;

int BspNode::mailID = 1;

/** Evaluates split plane classification with respect to the plane's 
	contribution for a minimum number splits in the tree.
*/
float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
/** Evaluates split plane classification with respect to the plane's 
	contribution for a balanced tree.
*/
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.
*/
float BspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
/** Evaluates split plane classification with respect to the plane's 
	contribution for balanced rays.
*/
float BspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};

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

BspNode::BspNode(): 
mParent(NULL)
{}

BspNode::BspNode(BspInterior *parent): 
mParent(parent)
{}


bool BspNode::IsRoot() const 
{
	return mParent == NULL;
}

BspInterior *BspNode::GetParent() 
{ 
	return mParent; 
}

void BspNode::SetParent(BspInterior *parent)
{
	mParent = parent;
}

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


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

bool BspInterior::IsLeaf() const
{ 
	return false; 
}

BspNode *BspInterior::GetBack() 
{
	return mBack;
}

BspNode *BspInterior::GetFront() 
{
	return mFront;
}

Plane3 *BspInterior::GetPlane()
{
	return &mPlane;
}

void BspInterior::ReplaceChildLink(BspNode *oldChild, BspNode *newChild) 
{
	if (mBack == oldChild)
		mBack = newChild;
	else
		mFront = newChild;
}

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

int BspInterior::SplitPolygons(PolygonContainer &polys, 
							   PolygonContainer &frontPolys, 
							   PolygonContainer &backPolys, 
							   PolygonContainer &coincident)
{
	Polygon3 *splitPoly = NULL;

	int splits = 0;

#ifdef _Debug
	Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
#endif
	while (!polys.empty())
	{
		Polygon3 *poly = polys.back();
		polys.pop_back();

		//Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";

		// classify polygon
		const int cf = poly->ClassifyPlane(mPlane);

		Polygon3 *front_piece = NULL;
		Polygon3 *back_piece = NULL;
	
		VertexContainer splitVertices;

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

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

				//-- inherit rays from parent polygon for blocked ray criterium
				poly->InheritRays(*front_piece, *back_piece);
				//Debug << "p: " << poly->mPiercingRays.size() << " f: " << front_piece->mPiercingRays.size() << " b: " << back_piece->mPiercingRays.size() << endl;
				
				// check if polygons still valid 
				if (front_piece->Valid())
					frontPolys.push_back(front_piece);
				else
					DEL_PTR(front_piece);
				
				if (back_piece->Valid())
					backPolys.push_back(back_piece);
				else				
					DEL_PTR(back_piece);
				
#ifdef _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;
}

/****************************************************************/
/*                  class BspLeaf implementation                */
/****************************************************************/
BspLeaf::BspLeaf(): mViewCell(NULL)
{
}

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

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

BspLeaf::BspLeaf(BspInterior *parent, BspViewCell *viewCell): 
BspNode(parent), mViewCell(viewCell)
{
}

BspViewCell *BspLeaf::GetViewCell() const
{
	return mViewCell;
}

void BspLeaf::SetViewCell(BspViewCell *viewCell)
{
	mViewCell = viewCell;
}

bool BspLeaf::IsLeaf() const 
{ 
	return true; 
}

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

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

	// add contributions from samples to the PVS
	for (it = rays.begin(); it != it_end; ++ it)
	{
		int contribution = 0;
		Ray *ray = (*it)->mRay;
			
		if (!ray->intersections.empty())
		{	
            //contribution += mViewCell->GetPvs().AddSample(ray->intersections[0].mObject);
		}

		//contribution += mViewCell->GetPvs().AddSample(ray->sourceObject.mObject);

		if (contribution > 0)
		{
			sampleContributions += contribution;
			++ contributingSamples;
		}

		// warning: not ordered
		ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
	}
}


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

BspTree::BspTree(BspViewCell *viewCell): 
mRootCell(viewCell), mRoot(NULL), mGenerateViewCells(false),
mStorePiercingRays(true)
{
	Randomize(); // initialise random generator for heuristics
}
 
const BspTreeStatistics &BspTree::GetStatistics() const 
{
	return mStat;
}

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

	app << setprecision(4);

	app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;

	app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;

	app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;

	app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;

	app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;

	app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;

	app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgBspLeaves() << endl;

	app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxBspLeaves << endl;
	
	app << "===== END OF BspViewCells statistics ==========\n";
}

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

	app << setprecision(4);

	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_SPLITS ( Number of splits )\n" << splits << "\n";

	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
	maxDepthNodes * 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_INPUT_POLYGONS (number of input polygons )\n" <<	polys << endl;

	//app << "#N_PVS: " << pvs << endl;

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


BspTree::~BspTree()
{
	std::stack<BspNode *> tStack;

	tStack.push(mRoot);

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

	    tStack.pop();
	
		if (!node->IsLeaf())
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

			// push the children on the stack (there are always two children)
			tStack.push(interior->GetBack());
			tStack.push(interior->GetFront());
		}

		DEL_PTR(node);
	}
}

void BspTree::InsertViewCell(ViewCell *viewCell)
{
	PolygonContainer *polys = new PolygonContainer();

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

	InsertPolygons(polys);
}

void BspTree::InsertPolygons(PolygonContainer *polys)
{	
	std::stack<BspTraversalData> tStack;
	
	// traverse existing tree or create new tree
    if (!mRoot)
		mRoot = new BspLeaf();

	tStack.push(BspTraversalData(mRoot, polys, 0, mRootCell, new BoundedRayContainer(), 0, 
								 mBox.SurfaceArea(), 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 += interior->SplitPolygons(*tData.mPolygons,
												  *frontPolys, 
												  *backPolys,
												  coincident);
				
				// extract view cells associated with the split polygons
				ViewCell *frontViewCell = mRootCell;
				ViewCell *backViewCell = mRootCell;
			
				BspTraversalData frontData(interior->GetFront(), 
										   frontPolys, 
										   tData.mDepth + 1,
										   mRootCell,	
										   tData.mRays,
										   tData.mPvs,
										   mBox.SurfaceArea(),
										   new BspNodeGeometry());

				BspTraversalData backData(interior->GetBack(), 
										  backPolys,
										  tData.mDepth + 1,
										  mRootCell,	
										  tData.mRays,
										  tData.mPvs,
										  mBox.SurfaceArea(),
										  new BspNodeGeometry());

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

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

				mStat.splits += splits;

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

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

int BspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
{
	FaceContainer::const_iterator fi;
	
	// copy the face data to polygons
	for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
	{
		Polygon3 *poly = new Polygon3((*fi), mesh);
		
		if (poly->Valid())
		{
			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)
	{
		if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
		{
			mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
			polysSize += AddMeshToPolygons(viewCells[i]->GetMesh(), polys, viewCells[i]);
		}
	}

	return polysSize;
}

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

		switch (object->Type()) // extract the meshes
		{
		case Intersectable::MESH_INSTANCE:
			mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
			break;
		case Intersectable::VIEW_CELL:
			mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
			break;
			// TODO: handle transformed mesh instances
		default:
			Debug << "intersectable type not supported" << endl;
			break;
		}
		
        if (mesh) // copy the mesh data to polygons
		{
			mBox.Include(object->GetBox()); // add to BSP tree aabb
			AddMeshToPolygons(mesh, polys, mRootCell);
		}
	}

	return (int)polys.size();
}

void BspTree::Construct(const ViewCellContainer &viewCells)
{
	mStat.nodes = 1;
	mBox.Initialize(); 	// initialise bsp tree bounding box

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

	// construct tree from the view cell polygons
	Construct(polys, new BoundedRayContainer());
}


void BspTree::Construct(const ObjectContainer &objects)
{
	mStat.nodes = 1;
	mBox.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::Construct(const RayContainer &sampleRays)
{
	mStat.nodes = 1;
	mBox.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;

	//-- extract polygons intersected by the rays
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		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);
			Face *face = obj->GetMesh()->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 (sSplitPlaneStrategy & BLOCKED_RAYS)
					(*it).second->mPiercingRays.push_back(ray);
			} 
			else
			{	//store rays if needed for heuristics
				Polygon3 *poly = new Polygon3(face, obj->GetMesh());
				poly->mParent = obj;
				polys->push_back(poly);

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

				facePolyMap[face] = poly;
			}
		}
	}
	
	facePolyMap.clear();

	// compue bounding box
	Polygon3::IncludeInBox(*polys, mBox);

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

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

	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(PolygonContainer *polys, BoundedRayContainer *rays)
{
	std::stack<BspTraversalData> tStack;

	mRoot = new BspLeaf();

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

	BspTraversalData tData(mRoot, polys, 0, mRootCell, rays, 
						   ComputePvsSize(*rays), cell->GetArea(), cell);

	tStack.push(tData);

	long startTime = GetTime();
	cout << "**** Contructing bsp tree ****\n";

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

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

	cout << "**** Finished tree construction ****\n";
	Debug << "BSP tree contruction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
}


BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
{
	//-- terminate traversal  
	if (((int)tData.mPolygons->size() <= sTermMaxPolygons) || 
		((int)tData.mRays->size() <= sTermMaxRays) ||
		(tData.mPvs <= sTermMinPvs) ||
		 (tData.mDepth >= sTermMaxDepth))
		
	{
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
	
		BspViewCell *viewCell;

		// generate new view cell for each leaf
		if (mGenerateViewCells)
		{
			viewCell = dynamic_cast<BspViewCell *>(ViewCell::Generate());
		}
		else
		{
			// add view cell to leaf
			viewCell = dynamic_cast<BspViewCell *>(tData.mViewCell);
		}

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

		//-- add pvs
		if (viewCell != mRootCell)
		{
			int conSamp = 0, sampCon = 0;
			leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
			
			mStat.contributingSamples += conSamp;
			mStat.sampleContributions += sampCon;
		}

		EvaluateLeafStats(tData);
		
		//-- clean up
		
		// discard polygons
		CLEAR_CONTAINER(*tData.mPolygons);
		// discard rays
		CLEAR_CONTAINER(*tData.mRays);

		DEL_PTR(tData.mPolygons);
		DEL_PTR(tData.mRays);
		DEL_PTR(tData.mCell);
		
		return leaf;
	}

	//-- continue subdivision

	PolygonContainer coincident;
	
	PolygonContainer *frontPolys = new PolygonContainer();
	PolygonContainer *backPolys = new PolygonContainer();

	BoundedRayContainer *frontRays = new BoundedRayContainer();
	BoundedRayContainer *backRays = new BoundedRayContainer();
	
	BspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
								new BoundedRayContainer(), 0, NULL, 0);
	BspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
							   new BoundedRayContainer(), 0, NULL, 0);

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

#ifdef _DEBUG	
	if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
	{	for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
			Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
		Debug << endl;}
#endif

	// extract view cells from coincident polygons according to plane normal
    // only if front or back polygons are empty
	if (!mGenerateViewCells)
	{
   		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.mCell);		

	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<ViewCell *>((*it)->mParent);
			foundBack = true;
		}
		else
		{
			frontData.mViewCell = dynamic_cast<ViewCell *>((*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, frontData, backData)); 

#ifdef _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.splits += interior->SplitPolygons(*tData.mPolygons, 
		                                    *frontData.mPolygons, 
											*backData.mPolygons, 
											coincident);

	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->mFront;
	backData.mNode = interior->mBack;
	
	return interior;
}

void BspTree::SortSplitCandidates(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;

	SortSplitCandidates(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 + sSplitBorder * (maxBox - minBox);
	float maxBand = minBox + (1.0f - sSplitBorder) * (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);

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

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


#if 0
  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
	Polygon3::IncludeInBox(polys, box);
        
	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 >= sMaxCostRatio)
		return false;

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

	return true;
}

Plane3 BspTree::SelectPlane(BspLeaf *leaf, 
							BspTraversalData &data,
							BspTraversalData &frontData,
							BspTraversalData &backData)
{
	if (data.mPolygons->empty() && data.mRays->empty())
	{
		Debug << "Warning: No autopartition polygon candidate available\n";
	
		// return axis aligned split
		AxisAlignedBox3 box;
		box.Initialize();
	
		// create bounding box of region
		Polygon3::IncludeInBox(*data.mPolygons, 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;
		return Plane3(norm, position);
	}
	
	if ((sSplitPlaneStrategy & AXIS_ALIGNED) &&
		((int)data.mPolygons->size() > sTermMaxPolysForAxisAligned) &&
		((int)data.mRays->size() > sTermMaxRaysForAxisAligned) &&
		((sTermMaxObjectsForAxisAligned < 0) || 
		  (Polygon3::ParentObjectsSize(*data.mPolygons) > sTermMaxObjectsForAxisAligned)))
	{
		Plane3 plane;
		if (SelectAxisAlignedPlane(plane, *data.mPolygons))
			return plane;
	}

	// simplest strategy: just take next polygon
	if (sSplitPlaneStrategy & RANDOM_POLYGON)
	{
		Polygon3 *nextPoly = (*data.mPolygons)[Random((int)data.mPolygons->size())];
		return nextPoly->GetSupportingPlane();
	}

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

Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf,
									  BspTraversalData &data,
									  BspTraversalData &frontData,
									  BspTraversalData &backData)
{
	float lowestCost = MAX_FLOAT;
	Plane3 bestPlane;
	
	int limit = sMaxPolyCandidates > 0 ? 
		Min((int)data.mPolygons->size(), sMaxPolyCandidates) : (int)data.mPolygons->size();
	
	int candidateIdx = limit;

	for (int i = 0; i < limit; ++ i)
	{
		candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
		
		Polygon3 *poly = (*data.mPolygons)[candidateIdx];
		// evaluate current candidate
		const float candidateCost = 
			SplitPlaneCost(poly->GetSupportingPlane(), data, frontData, backData);

		Debug << "cost: " << candidateCost << ", lowest: " << lowestCost << endl;

		if (candidateCost < lowestCost)
		{
			bestPlane = poly->GetSupportingPlane();
			lowestCost = candidateCost;
		}
	}
	
	//limit = maxRaysCandidates > 0 ? Min((int)rays.size(), maxRayCandidates) : (int)rays.size();	
	const BoundedRayContainer *rays = data.mRays;

	for (int i = 0; i < sMaxRayCandidates; ++ i)
	{
		Plane3 plane;

		if (1)
		{
			Vector3 pt[3];
			int idx[3];
			
			for (int j = 0; j < 3; j ++)
			{
				idx[j] = Random((int)rays->size() * 2);
				
				BoundedRay *bRay = (*rays)[idx[j]];
				Ray *ray = bRay->mRay;

				if (idx[j] < (int)rays->size())
					pt[j] = ray->Extrap(bRay->mMinT);
				else
				{
					idx[j] -= (int)rays->size();
					pt[j] = ray->Extrap(bRay->mMaxT);
				}
			}	
		
			plane = Plane3(pt[0], pt[1], pt[2]);
		}
		else
		{
			candidateIdx = Random((int)rays->size());
			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();
			
			plane = Plane3(normal, pt);
		}
		
		const float candidateCost = 
			SplitPlaneCost(plane, data, frontData, backData);

		if (candidateCost < lowestCost)
		{
			Debug << "choose ray plane: " << lowestCost << endl;
			bestPlane = plane;
			
			lowestCost = candidateCost;
		}
		else
			Debug << "ray cost: " << candidateCost << ", lowest: " << lowestCost << endl;
	}

	//Debug << "Plane lowest cost: " << lowestCost << endl;
	return bestPlane;
}

int BspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
{
	const int candidateIdx = Random(currentIdx --);

	// swap candidates to avoid testing same plane 2 times
	std::swap(polys[currentIdx], polys[candidateIdx]);
	
	return currentIdx;
	//return Random((int)polys.size());
}

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;

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

	for (it = polys.begin(); it != it_end; ++ it)
	{
		const int classification = (*it)->ClassifyPlane(candidatePlane);

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

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

		// assign view cells to back or front according to classificaion
		if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
		{
			MeshInstance *viewCell = (*it)->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;
				}
			}
		}
	}

	// all values should be approx. between 0 and 1 so they can be combined
	// and scaled with the factors according to their importance
	if (sSplitPlaneStrategy & BALANCED_POLYS)
		val += sBalancedPolysFactor * fabs(sumBalancedPolys) / (float)polys.size();
	
	if (sSplitPlaneStrategy & LEAST_SPLITS)    
		val += sLeastSplitsFactor * sumSplits / (float)polys.size();

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

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

	if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
		if (totalViewCells != 0)
			val += sBalancedViewCellsFactor * fabs(sumBalancedViewCells) / (float)totalViewCells;
	
	return val;
}

bool BspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
{
	maxT = 1e6;
	minT = 0;
	
	// test with tree bounding box
	if (!mBox.GetMinMaxT(ray, &minT, &maxT))
		return false;

	if (minT < 0) // start ray from origin
		minT = 0;

	// bound ray or line segment
	if ((ray.GetType() == Ray::LOCAL_RAY) && 
	    !ray.intersections.empty() && 
		(ray.intersections[0].mT <= maxT))
	{
		maxT = ray.intersections[0].mT;
	}

	return true;
}

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

	float sumBalancedRays = 0;
	float sumRaySplits = 0;
	
	int backId = 0;
	int frontId = 0;
	int frontAndBackId = 0;

	frontData.mPvs = 0;
	backData.mPvs = 0;
	frontData.mArea = 0;
	backData.mArea = 0;

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

	if (sSplitPlaneStrategy & PVS)
	{
		// create three unique ids for pvs heuristics
		Intersectable::NewMail(); backId = ViewCell::sMailId;
		Intersectable::NewMail(); frontId = ViewCell::sMailId;
		Intersectable::NewMail(); frontAndBackId = ViewCell::sMailId;

		if (1) // use front and back cell areas to approximate volume
		{	
			// construct child geometry with regard to the candidate split plane
			frontData.mCell = new BspNodeGeometry();
			backData.mCell = new BspNodeGeometry();

			cell.SplitGeometry(*frontData.mCell, *backData.mCell, *this, candidatePlane);
		
			pFront = frontData.mArea = frontData.mCell->GetArea();
			pBack = backData.mArea = backData.mCell->GetArea();

			pOverall = area;
		}
		else
		{
			pOverall = (float)rays.size();
		}
	}
			
	BoundedRayContainer::const_iterator rit, rit_end = rays.end();

	for (rit = rays.begin(); rit != rays.end(); ++ rit)
	{
		Ray *ray = (*rit)->mRay;
		const float minT = (*rit)->mMinT;
		const float maxT = (*rit)->mMaxT;

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

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

		if (sSplitPlaneStrategy & PVS)
		{
			if (!ray->intersections.empty())
			{
				// in case the ray intersects an objcrs
				// assure that we only count a object 
				// once for the front and once for the back side of the plane
				AddToPvs(*ray->intersections[0].mObject, frontData.mPvs, backData.mPvs, 
						 cf, frontId, backId, frontAndBackId);
			}

			// the source object in the origin of the ray
			AddToPvs(*ray->sourceObject.mObject, frontData.mPvs, backData.mPvs, 
					 cf, frontId, backId, frontAndBackId);

			// use #rays to approximate volume
			if (0)
			{
				if ((cf == Ray::BACK) || (cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT))
					++ pFront;
				if ((cf == Ray::FRONT) || (cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT))
					++ pBack;
			}
		}
	}

	if ((sSplitPlaneStrategy & LEAST_RAY_SPLITS) && !rays.empty())
			val += sLeastRaySplitsFactor * sumRaySplits / (float)rays.size();

	if ((sSplitPlaneStrategy & BALANCED_RAYS) && !rays.empty())
			val += sBalancedRaysFactor * fabs(sumBalancedRays) / (float)rays.size();

	if ((sSplitPlaneStrategy & PVS) && area && pvs)
		val += sPvsFactor * (frontData.mPvs * pFront + (backData.mPvs * pBack)) /
			   (pOverall * (float)pvs * 2);

	Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
		  << " frontpvs: " << frontData.mPvs << " pFront: " << pFront 
		  << " backpvs: " << backData.mPvs << " pBack: " << pBack 
		  << " val " << val << endl;

	return val;
}

void BspTree::AddToPvs(Intersectable &obj,
					   int &frontPvs,
					   int &backPvs,
					   const int cf, 
					   const int frontId, 
					   const int backId, 
					   const int frontAndBackId) const
{
	if (cf == Ray::COINCIDENT) //TODO: really belongs to no pvs?
		return;

	if (cf == Ray::FRONT)
	{
		if ((obj.mMailbox != frontId) && 
			(obj.mMailbox != frontAndBackId))
		{
			++ frontPvs;

			if (obj.mMailbox != backId)
				obj.mMailbox = frontId;
			else
				obj.mMailbox = frontAndBackId;					
		}
	}
	else if (cf == Ray::BACK)
	{
		if ((obj.mMailbox != backId) &&
			(obj.mMailbox != frontAndBackId))
		{
			++ backPvs;

			if (obj.mMailbox != frontId)
				obj.mMailbox = backId;
			else
				obj.mMailbox = frontAndBackId; 
		}
	}
	// object belongs to both pvs
	else if ((cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT))
	{
		if (obj.mMailbox !=  frontAndBackId)
		{
			if (obj.mMailbox != frontId)
				++ frontPvs;
			if (obj.mMailbox != backId)
				++ backPvs;
		
			obj.mMailbox = frontAndBackId;
		}
	}
}

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

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

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

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

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

void BspTree::ParseEnvironment()
{
	//-- parse bsp cell tree construction method
	char constructionMethodStr[60];
	
	environment->GetStringValue("BspTree.Construction.input", constructionMethodStr);

	sConstructionMethod = FROM_INPUT_VIEW_CELLS;
	
	if (strcmp(constructionMethodStr, "fromViewCells") == 0)
		sConstructionMethod = FROM_INPUT_VIEW_CELLS;
	else if (strcmp(constructionMethodStr, "fromSceneGeometry") == 0)
		sConstructionMethod = FROM_SCENE_GEOMETRY;
	else if (strcmp(constructionMethodStr, "fromRays") == 0)
		sConstructionMethod = FROM_RAYS;
	else 
	{
		cerr << "Wrong construction method " << constructionMethodStr << endl;
		exit(1);
    }

	Debug << "Construction method: " << constructionMethodStr << endl;

	//-- termination criteria for autopartition
	environment->GetIntValue("BspTree.Termination.maxDepth", sTermMaxDepth);
	environment->GetIntValue("BspTree.Termination.minPvs", sTermMinPvs);
	environment->GetIntValue("BspTree.Termination.maxPolygons", sTermMaxPolygons);
	environment->GetIntValue("BspTree.Termination.maxRays", sTermMaxRays);

	//-- termination criteria for axis aligned split
	environment->GetFloatValue("BspTree.Termination.AxisAligned.ct_div_ci", sCt_div_ci);
	environment->GetFloatValue("BspTree.Termination.AxisAligned.maxCostRatio", sMaxCostRatio);
	environment->GetIntValue("BspTree.Termination.AxisAligned.maxPolys", 
							 sTermMaxPolysForAxisAligned);
	environment->GetIntValue("BspTree.Termination.AxisAligned.maxRays", 
							 sTermMaxRaysForAxisAligned);
	environment->GetIntValue("BspTree.Termination.AxisAligned.maxObjects", 
							 sTermMaxObjectsForAxisAligned);
	//-- partition criteria
	environment->GetIntValue("BspTree.maxPolyCandidates", sMaxPolyCandidates);
	environment->GetIntValue("BspTree.maxRayCandidates", sMaxRayCandidates);
	environment->GetIntValue("BspTree.splitPlaneStrategy", sSplitPlaneStrategy);
	environment->GetFloatValue("BspTree.AxisAligned.splitBorder", sSplitBorder);

	environment->GetFloatValue("BspTree.Construction.sideTolerance", Vector3::sDistTolerance);
	Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;

	// post processing stuff
	environment->GetIntValue("ViewCells.PostProcessing.minPvsDif", sMinPvsDif);
	environment->GetIntValue("ViewCells.PostProcessing.minPvs", sMinPvs);
	environment->GetIntValue("ViewCells.PostProcessing.maxPvs", sMaxPvs);

    Debug << "BSP max depth: " << sTermMaxDepth << endl;
	Debug << "BSP min PVS: " << sTermMinPvs << endl;
	Debug << "BSP max polys: " << sTermMaxPolygons << endl;
	Debug << "BSP max rays: " << sTermMaxRays << endl;
	Debug << "BSP max polygon candidates: " << sMaxPolyCandidates << endl;
	Debug << "BSP max plane candidates: " << sMaxRayCandidates << endl;

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

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

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

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

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

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

	// store minimal and maximal pvs
	/*if (data.mPvs > mStat.pvs)
		mStat.pvs = data.mPvs;
	
	if (data.mPvs < mStat.pvs)
		mStat.pvs = data.mPvs;*/

	// accumulate depth to compute average
	mStat.accumDepth += data.mDepth;
	
//#ifdef _DEBUG
	Debug << "BSP stats: "
		  << "Depth: " << data.mDepth << " (max: " << sTermMaxDepth << "), "
		   << "PVS: " << data.mPvs << " (max: " << sTermMinPvs << "), "
		  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << sTermMaxPolygons << "), "
		  << "#rays: " << (int)data.mRays->size() << " (max: " << sTermMaxRays << "), " 
		  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << endl;
//#endif
}

int BspTree::CastRay(Ray &ray)
{
	int hits = 0;
  
	stack<BspRayTraversalData> tStack;
  
	float maxt, mint;

	if (!BoundRay(ray, mint, maxt))
		return 0;

	Intersectable::NewMail();

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

			int entSide = splitPlane->Side(entp);
			int extSide = splitPlane->Side(extp);

			Vector3 intersection;

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

bool BspTree::Export(const string filename)
{
	Exporter *exporter = Exporter::GetExporter(filename);

	if (exporter) 
	{
		exporter->ExportBspTree(*this);
		return true;
	}	

	return false;
}

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()) 
		{
			ViewCell *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->mFront);
			nodeStack.push(interior->mBack);
		}
	}
}

void BspTree::EvaluateViewCellsStats(BspViewCellsStatistics &stat) const
{
	stat.Reset();

	stack<BspNode *> nodeStack;
	nodeStack.push(mRoot);

	ViewCell::NewMail();

	// exclude root cell
	mRootCell->Mail();

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

		if (node->IsLeaf()) 
		{
			++ stat.bspLeaves;

			BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->mViewCell;

			if (!viewCell->Mailed()) 
			{
				viewCell->Mail();
				
				++ stat.viewCells;
				int pvsSize = viewCell->GetPvs().GetSize();

                stat.pvs += pvsSize;

				if (pvsSize < 2)
					++ stat.emptyPvs;

				if (pvsSize > stat.maxPvs)
					stat.maxPvs = pvsSize;

				if (pvsSize < stat.minPvs)
					stat.minPvs = pvsSize;

				if ((int)viewCell->mBspLeaves.size() > stat.maxBspLeaves)
					stat.maxBspLeaves = (int)viewCell->mBspLeaves.size();		
			}
		}
		else 
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

			nodeStack.push(interior->mFront);
			nodeStack.push(interior->mBack);
		}
	}
}

bool BspTree::MergeViewCells(BspLeaf *front, BspLeaf *back) const
{
	BspViewCell *viewCell = 
		dynamic_cast<BspViewCell *>(ViewCell::Merge(*front->mViewCell, *back->mViewCell));
	
	if (!viewCell)
		return false;

	// change view cells of all leaves associated with the
	// previous view cells

	BspViewCell *fVc = front->mViewCell;
	BspViewCell *bVc = back->mViewCell;

	vector<BspLeaf *> fLeaves = fVc->mBspLeaves;
	vector<BspLeaf *> bLeaves = bVc->mBspLeaves;

	vector<BspLeaf *>::const_iterator it;
	
	for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
	{
		(*it)->SetViewCell(viewCell);
		viewCell->mBspLeaves.push_back(*it);
	}
	for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
	{
		(*it)->SetViewCell(viewCell);
		viewCell->mBspLeaves.push_back(*it);
	}
	
	DEL_PTR(fVc);
	DEL_PTR(bVc);

	return true;
}

bool BspTree::ShouldMerge(BspLeaf *front, BspLeaf *back) const
{
	ViewCell *fvc = front->mViewCell;
	ViewCell *bvc = back->mViewCell;

	if ((fvc == mRootCell) || (bvc == mRootCell) || (fvc == bvc))
		return false;

	int fdiff = fvc->GetPvs().Diff(bvc->GetPvs());

	if (fvc->GetPvs().GetSize() + fdiff < sMaxPvs)
	{
		if ((fvc->GetPvs().GetSize() < sMinPvs) ||	
			(bvc->GetPvs().GetSize() < sMinPvs) ||
			((fdiff < sMinPvsDif) && (bvc->GetPvs().Diff(fvc->GetPvs()) < sMinPvsDif)))
		{
			return true;
		}
	}
	
	return false;
}

void BspTree::SetGenerateViewCells(int generateViewCells)
{
	mGenerateViewCells = generateViewCells;
}

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

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;
		rays.pop_back();
		
		const int cf = 
			ray->ClassifyPlane(plane, bRay->mMinT, bRay->mMaxT);
		
		ray->SetId(cf);

		switch (cf)
		{
		case Ray::COINCIDENT:
			//frontRays.push_back(bRay);
			DEL_PTR(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 Vector3 extp = ray->Extrap(bRay->mMaxT);
				const float t = plane.FindT(ray->GetLoc(), extp);
				
				const float newT = t * bRay->mMaxT;

				frontRays.push_back(new BoundedRay(ray, bRay->mMinT, newT));
				backRays.push_back(new BoundedRay(ray, newT, bRay->mMaxT));

				DEL_PTR(bRay);

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

				backRays.push_back(new BoundedRay(ray, bRay->mMinT, newT));
				frontRays.push_back(new BoundedRay(ray, newT, bRay->mMaxT));
				
				DEL_PTR(bRay);

				++ splits;
			}
			break;
		default:
			Debug << "Should not come here" << 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->mFront != lastNode)
				halfSpace.ReverseOrientation();

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

void BspTree::ConstructGeometry(BspNode *n, BspNodeGeometry &cell) const
{
	PolygonContainer polys;
	ConstructGeometry(n, polys);
	cell.mPolys = polys;
}

void BspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &cell) const
{
	vector<BspLeaf *> leaves = vc->mBspLeaves;

	vector<BspLeaf *>::const_iterator it, it_end = leaves.end();

	for (it = leaves.begin(); it != it_end; ++ it)
		ConstructGeometry(*it, cell);
}


void BspTree::ConstructGeometry(BspNode *n, PolygonContainer &cell) const
{
	vector<Plane3> halfSpaces;
	ExtractHalfSpaces(n, halfSpaces);

	PolygonContainer candidatePolys;

	// bounded planes are added to the polygons (reverse polygons
	// as they have to be outfacing
	for (int i = 0; i < (int)halfSpaces.size(); ++ i)
	{
		Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
		
		if (p->Valid())
		{
			candidatePolys.push_back(p->CreateReversePolygon());
			DEL_PTR(p);
		}
	}

	// add faces of bounding box (also could be faces of the cell)
	for (int i = 0; i < 6; ++ i)
	{
		VertexContainer vertices;
	
		for (int j = 0; j < 4; ++ j)
			vertices.push_back(mBox.GetFace(i).mVertices[j]);

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

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

			VertexContainer splitPts;
			Polygon3 *frontPoly, *backPoly;

			const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
			
			switch (cf)
			{
				case Polygon3::SPLIT:
					frontPoly = new Polygon3();
					backPoly = new Polygon3();

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

					DEL_PTR(candidatePolys[i]);

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

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


int BspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors, 
						   const bool onlyUnmailed) const
{
	PolygonContainer cell;

	ConstructGeometry(n, cell);

	stack<BspNode *> nodeStack;
	nodeStack.push(mRoot);
		
	// planes needed to verify that we found neighbor leaf.
	vector<Plane3> halfSpaces;
	ExtractHalfSpaces(n, halfSpaces);

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

		if (node->IsLeaf())
		{
            if (node != n && (!onlyUnmailed || !node->Mailed()))
			{
				// test all planes of current node if candidate really
				// is neighbour
				PolygonContainer neighborCandidate;
				ConstructGeometry(node, neighborCandidate);
				
				bool isAdjacent = true;
				for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
				{
					const int cf = 
						Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);

					if (cf == Polygon3::BACK_SIDE)
						isAdjacent = false;
				}

				if (isAdjacent)
					neighbors.push_back(dynamic_cast<BspLeaf *>(node));

				CLEAR_CONTAINER(neighborCandidate);
			}
		}
		else 
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);
	
			const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);

			if (cf == Polygon3::FRONT_SIDE)
				nodeStack.push(interior->mFront);
			else
				if (cf == Polygon3::BACK_SIDE)
					nodeStack.push(interior->mBack);
				else 
				{
					// random decision
					nodeStack.push(interior->mBack);
					nodeStack.push(interior->mFront);
				}
		}
	}
	
	CLEAR_CONTAINER(cell);
	return (int)neighbors.size();
}

BspLeaf *BspTree::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;
	
			PolygonContainer cell;

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

			const int cf = Polygon3::ClassifyPlane(cell, halfspace);

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

			nodeStack.push(next);
		}
	}
	
	return NULL;
}

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->mBack);
			else
				nodeStack.push(interior->mFront);

			mask = mask >> 1;
		}
	}
	
	return NULL;
}

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->Mailed())
		{
			ray->sourceObject.mObject->Mail();
			++ pvsSize;
		}
	}

	return pvsSize;
}

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

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

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

void BspNodeGeometry::SplitGeometry(BspNodeGeometry &front,
									BspNodeGeometry &back,
									const BspTree &tree,						 
									const Plane3 &splitPlane) const
{	
	// get cross section of new polygon
	Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);

	planePoly = SplitPolygon(planePoly, tree);

	//-- plane poly splits all other cell polygons
	for (int i = 0; i < (int)mPolys.size(); ++ i)
	{
		const int cf = mPolys[i]->ClassifyPlane(splitPlane);
			
		// split new polygon with all previous planes
		switch (cf)
		{
			case Polygon3::SPLIT:
				{
					Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);

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

					DEL_PTR(poly);

					if (frontPoly->Valid())
						front.mPolys.push_back(frontPoly);
					if (backPoly->Valid())
						back.mPolys.push_back(backPoly);
				}
				
				break;
			case Polygon3::BACK_SIDE:
				back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));			
				break;
			case Polygon3::FRONT_SIDE:
				front.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));	
				break;
			case Polygon3::COINCIDENT:
				//front.mPolys.push_back(CreateReversePolygon(mPolys[i]));
				//back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));
				break;
			default:
				break;
		}
	}

	//-- finally add the new polygon to the child cells
	if (planePoly)
	{
		// add polygon with normal pointing into positive half space to back cell
		back.mPolys.push_back(planePoly);
		// add polygon with reverse orientation to front cell
		front.mPolys.push_back(planePoly->CreateReversePolygon());
	}

	//Debug << "returning new geometry " << mPolys.size() << " CHILD: " << childCell->mPolys.size() << endl;
	//Debug << "old area " << this->GetArea() << " new: " << childCell->GetArea() << endl;
}

Polygon3 *BspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
										const BspTree &tree) const
{
	// polygon is split by all other planes
	for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
	{
		Plane3 plane = mPolys[i]->GetSupportingPlane();

		const int cf = 
			planePoly->ClassifyPlane(plane);
			
		// split new polygon with all previous planes
		switch (cf)
		{
			case Polygon3::SPLIT:
				{
					VertexContainer splitPts;
				
					Polygon3 *frontPoly = new Polygon3();
					Polygon3 *backPoly = new Polygon3();

					planePoly->Split(plane, *frontPoly, *backPoly, splitPts);
					DEL_PTR(planePoly);

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