#include "Plane3.h"
#include "VspBspTree.h"
#include "Mesh.h"
#include "common.h"
#include "ViewCell.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include <stack>
#include <time.h>
#include <iomanip>
#include "Exporter.h"
#include "Plane3.h"
#include "ViewCellBsp.h"

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


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


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

VspBspTree::VspBspTree(): 
mRoot(NULL),
mPvsUseArea(true)
{
	mRootCell = new BspViewCell();

	Randomize(); // initialise random generator for heuristics

	//-- termination criteria for autopartition
	environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
	environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
	environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
	environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);	
	environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
	environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);

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

	//-- termination criteria for axis aligned split
	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAaCtDivCi);
	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
	
	environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays", 
							 mTermMinRaysForAxisAligned);
	
	//-- partition criteria
	environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
	environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
	environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
	environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mSplitBorder);

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

    Debug << "BSP max depth: " << mTermMaxDepth << endl;
	Debug << "BSP min PVS: " << mTermMinPvs << endl;
	Debug << "BSP min area: " << mTermMinArea << endl;
	Debug << "BSP min 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_RAY_SPLITS)
		Debug << "least ray splits ";
	if (mSplitPlaneStrategy & BALANCED_RAYS)
		Debug << "balanced rays ";
	if (mSplitPlaneStrategy & PVS)
		Debug << "pvs";

	Debug << endl;
}


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


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

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

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

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

	return polysSize;
}

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

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

	return (int)polys.size();
}

void VspBspTree::Construct(const VssRayContainer &sampleRays)
{
    mStat.nodes = 1;
	mBox.Initialize(); 	// initialise BSP tree bounding box
	
	PolygonContainer polys;
	RayInfoContainer *rays = new RayInfoContainer();

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

	long startTime = GetTime();

	Debug << "**** Extracting polygons from rays ****\n";

	Intersectable::NewMail();

	//-- extract polygons intersected by the rays
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = *rit;
	
		if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
		{
			ray->mTerminationObject->Mail();
			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
			AddMeshToPolygons(obj->GetMesh(), polys, obj);
		}

		if (ray->mOriginObject && !ray->mOriginObject->Mailed())
		{
			ray->mOriginObject->Mail();
			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
			AddMeshToPolygons(obj->GetMesh(), polys, obj);
		}
	}

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

	//-- store rays
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = *rit;
		
		float minT, maxT;

		// TODO: not very efficient to implictly cast between rays types ...
		if (mBox.GetRaySegment(*ray, minT, maxT))
		{
			float len = ray->Length();
			
			if (!len) 
				len = Limits::Small;
			
			rays->push_back(RayInfo(ray, minT / len, maxT / len));
		}
	}

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

	// clean up polygons
	CLEAR_CONTAINER(polys);
}

void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
{
	std::stack<VspBspTraversalData> tStack;

	mRoot = new BspLeaf();

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

	VspBspTraversalData tData(mRoot, 
							  new PolygonContainer(polys), 
							  0,
							  rays, 
                              ComputePvsSize(*rays), 
							  geom->GetArea(), 
							  geom);

	tStack.push(tData);

	mStat.Start();
	cout << "Contructing vsp bsp tree ... ";

	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 << "s" << endl;
	}

	cout << "finished\n";

	mStat.Stop();
}

bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
{
	return 
		(((int)data.mRays->size() <= mTermMinRays) ||
		 (data.mPvs <= mTermMinPvs) ||
		 (data.mArea <= mTermMinArea) ||
		// (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
		 (data.mDepth >= mTermMaxDepth));
}

BspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, 
							   VspBspTraversalData &tData)
{
	//-- terminate traversal  
	if (TerminationCriteriaMet(tData))		
	{
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
	
		BspViewCell *viewCell = new BspViewCell();
		
		leaf->SetViewCell(viewCell);
		viewCell->mBspLeaves.push_back(leaf);

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

		EvaluateLeafStats(tData);
	
		//-- clean up
		
		DEL_PTR(tData.mPolygons);
		DEL_PTR(tData.mRays);
		DEL_PTR(tData.mGeometry);

		return leaf;
	}

	//-- continue subdivision
	PolygonContainer coincident;
	
	VspBspTraversalData tFrontData(new PolygonContainer(), 
								   tData.mDepth + 1,
								   new RayInfoContainer(), 
								   new BspNodeGeometry());

	VspBspTraversalData tBackData(new PolygonContainer(), 
						 		  tData.mDepth + 1,
								  new RayInfoContainer(), 
								  new BspNodeGeometry());

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

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

BspInterior *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
									   VspBspTraversalData &frontData,
									   VspBspTraversalData &backData,
									   PolygonContainer &coincident)
{
	mStat.nodes += 2;
	
	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
	
	// select subdivision plane
	BspInterior *interior = new BspInterior(SelectPlane(leaf, tData)); 

#ifdef _DEBUG
	Debug << interior << endl;
#endif
	
	// subdivide rays into front and back rays
	SplitRays(interior->GetPlane(), 
			  *tData.mRays, 
			  *frontData.mRays, 
			  *backData.mRays);
	
	// subdivide polygons with plane
	mStat.splits += 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->GetPlane(),
									   mBox,
									   mEpsilon);
	
		
		frontData.mArea = frontData.mGeometry->GetArea();
		backData.mArea = backData.mGeometry->GetArea();
	}

	// compute accumulated ray length
	//frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
	//backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);

	//-- 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();
	
	//DEL_PTR(leaf);
	return interior;
}

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

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

	ViewCell *vc = leaf->GetViewCell();

	// add contributions from samples to the PVS
	for (it = rays.begin(); it != it_end; ++ it)
	{
		int contribution = 0;
		VssRay *ray = (*it).mRay;
			
		if (ray->mTerminationObject)
			contribution += vc->GetPvs().AddSample(ray->mTerminationObject);
		
		if (ray->mOriginObject)
			contribution += vc->GetPvs().AddSample(ray->mOriginObject);

		if (contribution)
		{
			sampleContributions += contribution;
			++ contributingSamples;
		}
		//leaf->mVssRays.push_back(ray);
	}
}

void VspBspTree::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 VspBspTree::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 + 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);

			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 = mAaCtDivCi + 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 VspBspTree::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 >= mMaxCostRatio)
		return false;

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

	return true;
}

Plane3 VspBspTree::SelectPlane(BspLeaf *leaf, VspBspTraversalData &data)
{
	if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
		((int)data.mRays->size() > mTermMinRaysForAxisAligned))
	{
		Plane3 plane;
		if (SelectAxisAlignedPlane(plane, *data.mPolygons)) // TODO: should be rays
			return plane;
	}

	// simplest strategy: just take next polygon
	if (mSplitPlaneStrategy & RANDOM_POLYGON)
	{
        if (!data.mPolygons->empty())
		{
			const int randIdx = Random((int)data.mPolygons->size());
			Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
			return nextPoly->GetSupportingPlane();
		}
		else
		{
			//-- choose plane on midpoint of a ray
			const int candidateIdx = Random((int)data.mRays->size());
									
			const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
			const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();

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

			const Vector3 normal = (*data.mRays)[candidateIdx].mRay->GetDir();
			
			return Plane3(normal, pt);
		}

		return Plane3();
	}

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

Plane3 VspBspTree::SelectPlaneHeuristics(BspLeaf *leaf, VspBspTraversalData &data)
{
	float lowestCost = MAX_FLOAT;
	Plane3 bestPlane;
	Plane3 plane;

	int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
	
	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);

		if (candidateCost < lowestCost)
		{
			bestPlane = poly->GetSupportingPlane();
			lowestCost = candidateCost;
		}
	}
	
	//Debug << "lowest: " << lowestCost << endl;

	//-- choose candidate planes extracted from rays
	// we currently use different two methods chosen with
	// equal probability
	
	// take 3 ray endpoints, where two are minimum and one a maximum
	// point or the other way round
	for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
	{
		candidateIdx = Random((int)data.mRays->size());
	
		RayInfo rayInf = (*data.mRays)[candidateIdx];

		const Vector3 minPt = rayInf.ExtrapOrigin();
		const Vector3 maxPt = rayInf.ExtrapTermination();

		const Vector3 pt = (maxPt + minPt) * 0.5;
		const Vector3 normal = Normalize(rayInf.mRay->GetDir());

		plane = Plane3(normal, pt);

		const float candidateCost = SplitPlaneCost(plane, data);

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

	// take plane normal as plane normal and the midpoint of the ray.
	// PROBLEM: does not resemble any point where visibility is likely to change
	//Debug << "lowest2: " << lowestCost << endl;
	for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
	{
		Vector3 pt[3];
		int idx[3];
		int cmaxT = 0;
		int cminT = 0;
		bool chooseMin = false;

		for (int j = 0; j < 3; j ++)
		{
			idx[j] = Random((int)data.mRays->size() * 2);
				
			if (idx[j] >= (int)data.mRays->size())
			{
				idx[j] -= (int)data.mRays->size();
				
				chooseMin = (cminT < 2);
			}
			else
				chooseMin = (cmaxT < 2);

			RayInfo rayInf = (*data.mRays)[idx[j]];
			pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
		}	
			
		plane = Plane3(pt[0], pt[1], pt[2]);

		const float candidateCost = SplitPlaneCost(plane, data);

		if (candidateCost < lowestCost)
		{
			//Debug << "choose ray plane 2: " << candidateCost << endl;
			bestPlane = plane;
			
			lowestCost = candidateCost;
		}
	}	

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

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


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

float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
								 const VspBspTraversalData &data)
{
	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();
	
		if (mPvsUseArea) // use front and back cell areas to approximate volume
		{	
			// construct child geometry with regard to the candidate split plane
			BspNodeGeometry frontCell;
			BspNodeGeometry backCell;
		
			data.mGeometry->SplitGeometry(frontCell, 
										  backCell, 
										  candidatePlane,
										  mBox,
										  mEpsilon);
		
			pFront = frontCell.GetArea();
			pBack = backCell.GetArea();

			pOverall = data.mArea;
		}
	}
		
	int limit;
	bool useRand;

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

	for (int i = 0; i < limit; ++ i)
	{
		const int testIdx = useRand ? Random(limit) : i;
		RayInfo rayInf = (*data.mRays)[testIdx];

		VssRay *ray = rayInf.mRay;
		const int cf = rayInf.ComputeRayIntersection(candidatePlane, ray->mT);

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

		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
			AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
			
			// add the source object
			AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
			
			// use number or length of rays to approximate volume
			if (!mPvsUseArea)
			{
				float len = 1;

				if (pvsUseLen) // use length of rays
					len = rayInf.SqrSegmentLength();
			
				pOverall += len;

				if (cf == 1)
					pFront += len;
				if (cf == -1)
					pBack += len;
				if (cf == 0)
				{
					// use length of rays to approximate volume
					if (pvsUseLen) 
					{
						float newLen = len * 
							(rayInf.GetMaxT() - rayInf.mRay->mT) / 
							(rayInf.GetMaxT() - rayInf.GetMinT());
		
						if (candidatePlane.Side(rayInf.ExtrapOrigin()) <= 0)
						{
							pBack += newLen;
							pFront += len - newLen;
						}
						else
						{
							pFront += newLen;
							pBack += len - newLen;
						}
					}
					else
					{
						++ pFront;
						++ pBack;
					}
				}
			}
		}
	}

	const float raysSize = (float)data.mRays->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)data.mPvs * 2.0f + Limits::Small;

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

		// give penalty to unbalanced split
		if (0)
		if (((pFront * 0.2 + Limits::Small) > pBack) || 
			(pFront < (pBack * 0.2 + Limits::Small)))
			val += 0.5;
	}

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

void VspBspTree::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
	if (cf >= 0)
	{
		if ((obj->mMailbox != sFrontId) && 
			(obj->mMailbox != sFrontAndBackId))
		{
			++ frontPvs;

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

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

void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
{
	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 VspBspTree::GetBoundingBox() const
{
	return mBox;
}

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

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

	// store maximal and minimal depth
	if (data.mDepth > 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;
	

	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.mGeometry->GetArea() <= mTermMinArea) 
		++ mStat.minAreaNodes;

#ifdef _DEBUG
	Debug << "BSP stats: "
		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
		  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
		  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
		  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), " 
		  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
		  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
#endif
}

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

	if (!mBox.GetRaySegment(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 = 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->GetViewCell()->Mailed())
			{
				//ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
				leaf->GetViewCell()->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 VspBspTree::Export(const string filename)
{
	Exporter *exporter = Exporter::GetExporter(filename);

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

	return false;
}

void VspBspTree::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)->GetViewCell();

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

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

void VspBspTree::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)->GetViewCell();

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

                stat.pvs += pvsSize;

				if (pvsSize < 1)
					++ 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->GetFront());
			nodeStack.push(interior->GetBack());
		}
	}
}

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

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

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

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

	return len;
}

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

	while (!rays.empty())
	{
		RayInfo bRay = rays.back();
		VssRay *ray = bRay.mRay;

		rays.pop_back();
	
		// get classification and receive new t
		const int cf = bRay.ComputeRayIntersection(plane, ray->mT);

		switch (cf)
		{
		case -1:
			backRays.push_back(bRay);
			break;
		case 1:
			frontRays.push_back(bRay);
			break;
		case 0: 
			//-- split ray

			// if start point behind plane
			if (plane.Side(bRay.ExtrapOrigin()) <= 0)
			{	
				backRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
				frontRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
			}
			else
			{
				frontRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
				backRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
			}
			break;
		default:
			Debug << "Should not come here 4" << endl;
			break;
		}
	}

	return splits;
}

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

	do
	{
		lastNode = n;

		// want to get planes defining geometry of this node => don't take
		// split plane of node itself
		n = n->GetParent();
		
		if (n)
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(n);
			Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();

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

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

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

void VspBspTree::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(mEpsilon))
		{
			candidatePolys.push_back(p->CreateReversePolygon());
			DEL_PTR(p);
		}
	}

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

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

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

			VertexContainer splitPts;
			Polygon3 *frontPoly, *backPoly;

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

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

					DEL_PTR(candidatePolys[i]);

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

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

void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &vcGeom) 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, vcGeom);
}

int VspBspTree::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],
												mEpsilon);

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

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

BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
{
    stack<BspNode *> nodeStack;
  	nodeStack.push(mRoot);
	
	int mask = rand();
  
	while (!nodeStack.empty()) 
	{
		BspNode *node = nodeStack.top();
		nodeStack.pop();
	  
		if (node->IsLeaf()) 
		{
			return dynamic_cast<BspLeaf *>(node);
		} 
		else 
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);
			
			BspNode *next;
	
			PolygonContainer cell;

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

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

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

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

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

	int mask = rand();
	
	while (!nodeStack.empty()) 
	{
		BspNode *node = nodeStack.top();
		nodeStack.pop();
		
		if (node->IsLeaf()) 
		{
			if ( (!onlyUnmailed || !node->Mailed()) )
				return dynamic_cast<BspLeaf *>(node);
		}
		else 
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

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

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

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

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

	Intersectable::NewMail();

	for (rit = rays.begin(); rit != rays.end(); ++ rit)
	{
		VssRay *ray = (*rit).mRay;
		
		if (ray->mOriginObject)
		{
			if (!ray->mOriginObject->Mailed())
			{
				ray->mOriginObject->Mail();
				++ pvsSize;
			}
		}
		if (ray->mTerminationObject)
		{
			if (!ray->mTerminationObject->Mailed())
			{
				ray->mTerminationObject->Mail();
				++ pvsSize;
			}
		}
	}

	return pvsSize;
}

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

BspViewCell *VspBspTree::GetRootCell() const
{
	return mRootCell;
}

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

	while (!polys.empty())
	{
		Polygon3 *poly = polys.back();
		polys.pop_back();

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

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

	return splits;
}