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

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


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

float BspMergeCandidate::sOverallCost = Limits::Small;

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

VspBspTree::VspBspTree():
mRoot(NULL),
mPvsUseArea(true),
mCostNormalizer(Limits::Small),
mViewCellsManager(NULL),
mStoreRays(true),
mOnlyDrivingAxis(false)
{
	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);
	environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
	environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);

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

	//-- max cost ratio for early tree termination
	environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);

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

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

	// maximum and minimum number of view cells
	environment->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
	environment->GetIntValue("VspBspTree.PostProcess.minViewCells", mMergeMinViewCells);
	environment->GetFloatValue("VspBspTree.PostProcess.maxCostRatio", mMergeMaxCostRatio);


	//-- debug output
	Debug << "******* VSP BSP options ******** " << endl;
    Debug << "max depth: " << mTermMaxDepth << endl;
	Debug << "min PVS: " << mTermMinPvs << endl;
	Debug << "min area: " << mTermMinArea << endl;
	Debug << "min rays: " << mTermMinRays << endl;
	Debug << "max ray contri: " << mTermMaxRayContribution << endl;
	//Debug << "VSP BSP mininam accumulated ray lenght: ", mTermMinAccRayLength) << endl;
	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	Debug << "miss tolerance: " << mTermMissTolerance << endl;
	Debug << "max view cells: " << mMaxViewCells << endl;
	Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
	Debug << "max plane candidates: " << mMaxRayCandidates << endl;

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

	mSplitCandidates = new vector<SortableEntry>;

	Debug << endl;
}


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


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

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

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

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

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

	int polysSize = 0;

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

	return polysSize;
}

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

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

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

        if (mesh) // copy the mesh data to polygons
		{
			mBox.Include(object->GetBox()); // add to BSP tree aabb
			AddMeshToPolygons(mesh, polys, 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();

	cout << "Extracting polygons from rays ... ";

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

	mTermMinArea *= mBox.SurfaceArea(); // normalize
	mStat.polys = (int)polys.size();

	cout << "finished" << endl;

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

	Construct(polys, rays);

	// clean up polygons
	CLEAR_CONTAINER(polys);
}

void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
{
	VspBspTraversalStack 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 << "VSP BSP tree construction time spent at root: "
				  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl << 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) ||
		 (mStat.nodes / 2 + 1 >= mMaxViewCells) ||
		// (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
		 (data.mDepth >= mTermMaxDepth));
}

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

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

		VspBspTraversalData tFrontData;
		VspBspTraversalData tBackData;

		// create new interior node and two leaf nodes
		// or return leaf as it is (if maxCostRatio missed)
		newNode = SubdivideNode(tData, tFrontData, tBackData, coincident);

		if (!newNode->IsLeaf()) //-- continue subdivision
		{
			// push the children on the stack
			tStack.push(tFrontData);
			tStack.push(tBackData);

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

	//-- terminate traversal and create new view cell
	if (newNode->IsLeaf())
	{
		BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);

		// create new view cell for this leaf
		BspViewCell *viewCell = new BspViewCell();
		leaf->SetViewCell(viewCell);

		if (mStoreRays)
		{
			RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
			for (it = tData.mRays->begin(); it != it_end; ++ it)
				leaf->mVssRays.push_back((*it).mRay);
		}

		viewCell->mLeaves.push_back(leaf);
		viewCell->SetArea(tData.mArea);

		//-- update pvs
		int conSamp = 0, sampCon = 0;
		AddToPvs(leaf, *tData.mRays, conSamp, sampCon);

		mStat.contributingSamples += conSamp;
		mStat.sampleContributions += sampCon;

		EvaluateLeafStats(tData);
	}


	//-- cleanup
	tData.Clear();

	return newNode;
}

BspNode *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
								   VspBspTraversalData &frontData,
								   VspBspTraversalData &backData,
								   PolygonContainer &coincident)
{
	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);

	int maxCostMisses = tData.mMaxCostMisses;

	// select subdivision plane
	Plane3 splitPlane;
	if (!SelectPlane(splitPlane, leaf, tData))
	{
		++ maxCostMisses;

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

	mStat.nodes += 2;

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

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

	//-- the front and back traversal data is filled with the new values
	frontData.mPolygons = new PolygonContainer();
	frontData.mDepth = tData.mDepth + 1;
	frontData.mRays = new RayInfoContainer();
	frontData.mGeometry = new BspNodeGeometry();

	backData.mPolygons = new PolygonContainer();
	backData.mDepth = tData.mDepth + 1;
	backData.mRays = new RayInfoContainer();
	backData.mGeometry = new BspNodeGeometry();

	// subdivide rays
	SplitRays(interior->GetPlane(),
			  *tData.mRays,
			  *frontData.mRays,
			  *backData.mRays);

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


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

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

	// split 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 RayInfoContainer &rays, const int axis)
{
	mSplitCandidates->clear();

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

	mSplitCandidates->reserve(requestedSize);

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

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

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

float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
										  const AxisAlignedBox3 &box,
										  const int pvsSize,
										  const int &axis,
                                          float &position)
{
	int raysBack;
	int raysFront;
	int pvsBack;
	int pvsFront;

	SortSplitCandidates(rays, axis);

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

	int rl = 0, rr = (int)rays.size();
	int pl = 0, pr = pvsSize;

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

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

	float sum = rr*sizeBox;
	float minSum = 1e20f;

	Intersectable::NewMail();

	// set all object as belonging to the fron pvs
	for(RayInfoContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
	{
		if ((*ri).mRay->IsActive())
		{
			Intersectable *object = (*ri).mRay->mTerminationObject;

			if (object)
			{
				if (!object->Mailed())
				{
					object->Mail();
					object->mCounter = 1;
				}
				else
					++ object->mCounter;
			}
		}
	}

	Intersectable::NewMail();

	for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
		ci < mSplitCandidates->end(); ++ ci)
	{
		VssRay *ray;

		switch ((*ci).type)
		{
			case SortableEntry::ERayMin:
				{
					++ rl;
					ray = (VssRay *) (*ci).ray;

					Intersectable *object = ray->mTerminationObject;

					if (object && !object->Mailed())
					{
						object->Mail();
						++ pl;
					}
					break;
				}
			case SortableEntry::ERayMax:
				{
					-- rr;
					ray = (VssRay *) (*ci).ray;

					Intersectable *object = ray->mTerminationObject;

					if (object)
					{
						if (-- object->mCounter == 0)
							-- pr;
					}

					break;
				}
		}

		// Note: sufficient to compare size of bounding boxes of front and back side?
		if ((*ci).value > minBand && (*ci).value < maxBand)
		{
			sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);

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

			if (sum < minSum)
			{
				minSum = sum;
				position = (*ci).value;

				raysBack = rl;
				raysFront = rr;

				pvsBack = pl;
				pvsFront = pr;

			}
		}
	}

	float oldCost = (float)pvsSize;
	float newCost = mCtDivCi + minSum / sizeBox;
	float ratio = newCost / oldCost;

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

	return ratio;
}


float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
										 const VspBspTraversalData &tData)
{
	AxisAlignedBox3 box;
	box.Initialize();

	// create bounding box of region
	RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();

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

	const bool useCostHeuristics = false;

	//-- regular split
	float nPosition[3];
	float nCostRatio[3];
	int bestAxis = -1;

	const int sAxis = box.Size().DrivingAxis();

	for (int axis = 0; axis < 3; ++ axis)
	{
		if (!mOnlyDrivingAxis || axis == sAxis)
		{
			if (!useCostHeuristics)
			{
				nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
				Vector3 normal(0,0,0); normal[axis] = 1;

				nCostRatio[axis] = SplitPlaneCost(Plane3(normal, nPosition[axis]), tData);
			}
			else
			{
				nCostRatio[axis] =
					BestCostRatioHeuristics(*tData.mRays,
										    box,
											tData.mPvs,
											axis,
											nPosition[axis]);
			}

			if (bestAxis == -1)
			{
				bestAxis = axis;
			}

			else if (nCostRatio[axis] < nCostRatio[bestAxis])
			{
				/*Debug << "pvs front " << nPvsBack[axis]
					  << " pvs back " << nPvsFront[axis]
					  << " overall pvs " << leaf->GetPvsSize() << endl;*/
				bestAxis = axis;
			}

		}
	}

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

	return nCostRatio[bestAxis];
}


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

			plane = nextPoly->GetSupportingPlane();
			return true;
		}
		else
		{
			//-- choose plane on midpoint of a ray
			const int candidateIdx = (int)RandomValue(0, (Real)((int)data.mRays->size() - 1));

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

			plane = Plane3(normal, pt);
			return true;
		}

		return false;
	}

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


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

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

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

	return Plane3(normal, pt);
}


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	return Plane3(norm, planePt);
}

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

	const int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
	int maxIdx = (int)data.mPolygons->size();

	float candidateCost;

	for (int i = 0; i < limit; ++ i)
	{
		// assure that no index is taken twice
		const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
		//Debug << "current Idx: " << maxIdx << " cand idx " << candidateIdx << endl;

		Polygon3 *poly = (*data.mPolygons)[candidateIdx];

		// swap candidate to the end to avoid testing same plane
		std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);

		//Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];

		// evaluate current candidate
		candidateCost = SplitPlaneCost(poly->GetSupportingPlane(), data);

		if (candidateCost < lowestCost)
		{
			bestPlane = poly->GetSupportingPlane();
			lowestCost = candidateCost;
		}
	}

#if 0
	//-- choose candidate planes extracted from rays
	//-- different methods are available
	for (int i = 0; i < mMaxRayCandidates; ++ i)
	{
		plane = ChooseCandidatePlane3(*data.mRays);
		candidateCost = SplitPlaneCost(plane, data);

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

	// axis aligned splits
	candidateCost = SelectAxisAlignedPlane(plane, data);

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

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

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

	return true;
}


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) const
{
	float cost = 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 ? (int)RandomValue(0, (Real)(limit - 1)) : i;
		RayInfo rayInf = (*data.mRays)[testIdx];

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

		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() - t) /
							(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)
		cost += mLeastRaySplitsFactor * sumRaySplits / raysSize;

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

	// pvs criterium
	if (mSplitPlaneStrategy & PVS)
	{
		const float oldCost = pOverall * (float)data.mPvs + Limits::Small;
		cost += mPvsFactor * (frontPvs * pFront + backPvs * pBack) / oldCost;

		//cost += mPvsFactor * 0.5 * (frontPvs * pFront + backPvs * pBack) / oldCost;
		//Debug << "new cost: " << cost << " over" << frontPvs * pFront + backPvs * pBack << " old cost " << oldCost << endl;

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

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

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

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;

	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.mArea <= mTermMinArea)
	{
		//Debug << "area: " << data.mArea / mBox.SurfaceArea() << " min area: " << mTermMinArea / mBox.SurfaceArea() << endl;
		++ mStat.minAreaNodes;
	}
	// accumulate depth to compute average depth
	mStat.accumDepth += data.mDepth;

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


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();
		rays.pop_back();

		VssRay *ray = bRay.mRay;
		float t;

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

		switch (cf)
		{
		case -1:
			backRays.push_back(bRay);
			break;
		case 1:
			frontRays.push_back(bRay);
			break;
		case 0:
			//-- split ray
			//--  look if start point behind or in front of plane
			if (plane.Side(bRay.ExtrapOrigin()) <= 0)
			{
				backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
				frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
			}
			else
			{
				frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
				backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
			}
			break;
		default:
			Debug << "Should not come here 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->mLeaves;

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


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

	float mint = 0.0f, maxt = 1.0f;

	Intersectable::NewMail();

	Vector3 entp = origin;
	Vector3 extp = termination;

	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(origin, extp, &t);
			maxt *= t;

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

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

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

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


			BspRayTraversalData &s = tStack.top();

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

			tStack.pop();
		}
	}
	return hits;
}

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

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

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

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

	BspNode *p2 = n2;

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

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

	return 0; // never come here
}

BspNode *VspBspTree::CollapseTree(BspNode *node)
{
    if (node->IsLeaf())
		return node;

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

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

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

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

			BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);

			// replace a link from node's parent
			if (leaf->GetParent())
				leaf->GetParent()->ReplaceChildLink(node, leaf);

			delete interior;

			return leaf;
		}
	}

	return node;
}


void VspBspTree::RepairVcLeafLists()
{
	// list not valid anymore => clear
	stack<BspNode *> nodeStack;
	nodeStack.push(mRoot);

	ViewCell::NewMail();

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

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

			BspViewCell *viewCell = leaf->GetViewCell();

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

			viewCell->mLeaves.push_back(leaf);
		}
		else
		{
			BspInterior *interior = dynamic_cast<BspInterior *>(node);

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


int VspBspTree::MergeLeaves()
{
	vector<BspLeaf *> leaves;
	priority_queue<BspMergeCandidate> mergeQueue;

	// collect the leaves, e.g., the "voxels" that will build the view cells
	CollectLeaves(leaves);

	int vcSize = (int)leaves.size();
	int savedVcSize = vcSize;

	BspLeaf::NewMail();

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

	// find merge candidates and push them into queue
	for (it = leaves.begin(); it != it_end; ++ it)
	{
		BspLeaf *leaf = *it;

		// no leaf is part of two merge candidates
		if (!leaf->Mailed())
		{
			leaf->Mail();

			vector<BspLeaf *> neighbors;
			FindNeighbors(leaf, neighbors, true);

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

			for (nit = neighbors.begin(); nit != nit_end; ++ nit)
			{
				BspMergeCandidate mc = BspMergeCandidate(leaf, *nit);
				mergeQueue.push(mc);

				BspMergeCandidate::sOverallCost += mc.GetLeaf1Cost();
				BspMergeCandidate::sOverallCost += mc.GetLeaf2Cost();
			}
		}
	}

	int merged = 0;
	int mergedSiblings = 0;
	int accTreeDist = 0;
	int maxTreeDist = 0;
	const bool mergeStats = true;


	//-- use priority queue to merge leaf pairs
	while (!mergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
		   (mergeQueue.top().GetMergeCost() <
		    mMergeMaxCostRatio * BspMergeCandidate::sOverallCost))
	{
		//Debug << "mergecost: " << mergeQueue.top().GetMergeCost() / BspMergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
		BspMergeCandidate mc = mergeQueue.top();
		mergeQueue.pop();

		// both view cells equal!
		if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
			continue;

		if (mc.Valid())
		{
			ViewCell::NewMail();
			MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
			-- vcSize;
			// increase absolute merge cost
			BspMergeCandidate::sOverallCost += mc.GetMergeCost();

			//-- stats
			++ merged;

			if (mc.GetLeaf1()->IsSibling(mc.GetLeaf2()))
				++ mergedSiblings;

			if (mergeStats)
			{
				const int dist = TreeDistance(mc.GetLeaf1(), mc.GetLeaf2());

				if (dist > maxTreeDist)
					maxTreeDist = dist;

				accTreeDist += dist;
			}
		}
		// merge candidate not valid, because one of the leaves was already
		// merged with another one
		else
		{
			// validate and reinsert into queue
			mc.SetValid();
			mergeQueue.push(mc);
			//Debug << "validate " << mc.GetMergeCost() << endl;
		}
	}

	// collapse tree according to view cell partition
	CollapseTree(mRoot);
	// revalidate leaves
	RepairVcLeafLists();

	Debug << "Merged " << merged << " nodes of " << savedVcSize
		  << " (merged " << mergedSiblings << " siblings)" << endl;

	if (mergeStats)
	{
		Debug << "maximal tree distance: " << maxTreeDist << endl;
		Debug << "Avg tree distance: " << (float)accTreeDist / (float)merged << endl;
	}

	//TODO: should return sample contributions
	return merged;
}


bool VspBspTree::MergeViewCells(BspLeaf *l1, BspLeaf *l2)
{
	//-- change pointer to view cells of all leaves associated
	//-- with the previous view cells
	BspViewCell *fVc = l1->GetViewCell();
	BspViewCell *bVc = l2->GetViewCell();

	BspViewCell *vc = dynamic_cast<BspViewCell *>(
		mViewCellsManager->MergeViewCells(*fVc, *bVc));

	// if merge was unsuccessful
	if (!vc) return false;

	// set new size of view cell
	vc->SetArea(fVc->GetArea() + bVc->GetArea());

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

	vector<BspLeaf *>::const_iterator it;

	//-- change view cells of all the other leaves the view cell belongs to
	for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
	{
		(*it)->SetViewCell(vc);
		vc->mLeaves.push_back(*it);
	}

	for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
	{
		(*it)->SetViewCell(vc);
		vc->mLeaves.push_back(*it);
	}

	vc->Mail();

	// clean up old view cells
	DEL_PTR(fVc);
	DEL_PTR(bVc);

	return true;
}


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

/************************************************************************/
/*                BspMergeCandidate implementation                      */
/************************************************************************/


BspMergeCandidate::BspMergeCandidate(BspLeaf *l1, BspLeaf *l2):
mMergeCost(0),
mLeaf1(l1),
mLeaf2(l2),
mLeaf1Id(l1->GetViewCell()->mMailbox),
mLeaf2Id(l2->GetViewCell()->mMailbox)
{
	EvalMergeCost();
}

float BspMergeCandidate::GetLeaf1Cost() const
{
	BspViewCell *vc = mLeaf1->GetViewCell();
	return vc->GetPvs().GetSize() * vc->GetArea();
}

float BspMergeCandidate::GetLeaf2Cost() const
{
	BspViewCell *vc = mLeaf2->GetViewCell();
	return vc->GetPvs().GetSize() * vc->GetVolume();
}

void BspMergeCandidate::EvalMergeCost()
{
	//-- compute pvs difference
	BspViewCell *vc1 = mLeaf1->GetViewCell();
	BspViewCell *vc2 = mLeaf2->GetViewCell();

	const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
	const int vcPvs = diff1 + vc1->GetPvs().GetSize();

	//-- compute ratio of old cost
	//-- (added size of left and right view cell times pvs size)
	//-- to new rendering cost (size of merged view cell times pvs size)
	const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();

	const float newCost =
		(float)vcPvs * (vc1->GetArea() + vc2->GetArea());

	mMergeCost = newCost - oldCost;

//	if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
//		mMergeCost += 1.0;
}

void BspMergeCandidate::SetLeaf1(BspLeaf *l)
{
	mLeaf1 = l;
}

void BspMergeCandidate::SetLeaf2(BspLeaf *l)
{
	mLeaf2 = l;
}

BspLeaf *BspMergeCandidate::GetLeaf1()
{
	return mLeaf1;
}

BspLeaf *BspMergeCandidate::GetLeaf2()
{
	return mLeaf2;
}

bool BspMergeCandidate::Valid() const
{
	return
		(mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
		(mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
}

float BspMergeCandidate::GetMergeCost() const
{
	return mMergeCost;
}

void BspMergeCandidate::SetValid()
{
	mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
	mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;

	EvalMergeCost();
}