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

#include "ViewCell.h"
#include "Plane3.h"
#include "OspTree.h"
#include "Mesh.h"
#include "common.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include "Exporter.h"
#include "Plane3.h"
#include "ViewCellsManager.h"
#include "Beam.h"
#include "KdTree.h"
#include "KdIntersectable.h"
#include "VspTree.h"


namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0


//-- static members

OspTree *OspTree::OspSubdivisionCandidate::sOspTree = NULL;

// variable for debugging volume contribution for heuristics
static float debugVol;

static bool ViewCellHasMultipleReferences(Intersectable *obj, 
										  ViewCell *vc, 
										  bool checkOnlyMailed)
{
	MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);
//return false;
	if (vdata)
	{
		// more than one view cell sees this object inside different kd cells
		if (!checkOnlyMailed || !vdata->Mailed())
		{
			if (checkOnlyMailed)
				vdata->Mail();
			//Debug << "sumpdf: " << vdata->mSumPdf << endl;
			if (vdata->mSumPdf > 1.5f) 
				return true;
		}
	}

	return false;
}


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

	app << setprecision(4);

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

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

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

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

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

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

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

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

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

	//app << "#N_PMINOBJECTREFLEAVES  ( Percentage of leaves with minimal number of object ref)\n" 
		//<< minObjectNodes * 100 / (double)Leaves() << endl;

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

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

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

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

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

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

	app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" <<
		 maxObjectRefs << "\n";

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

	app << "========== END OF VspTree statistics ==========\n";
}


/******************************************************************/
/*                  class VspNode implementation                  */
/******************************************************************/


/*******************************************************************/
/*                  class OspTree implementation                   */
/*******************************************************************/


OspTree::OspTree():
mRoot(NULL),
mTimeStamp(1),
mCopyFromKdTree(false)
{
	ReadEnvironment();
	mSubdivisionCandidates = new vector<SortableEntry>;
}


OspTree::OspTree(const KdTree &kdTree)
{
	ReadEnvironment();

	// copy values from kd tree
	mCopyFromKdTree = true;
	mBoundingBox = kdTree.GetBox();
	mRoot = kdTree.GetRoot();

	mSubdivisionCandidates = new vector<SortableEntry>;
}


OspTree::~OspTree()
{
	// delete kd intersectables
	KdIntersectableMap::iterator it, it_end = mKdIntersectables.end();

	for (it = mKdIntersectables.begin(); it != mKdIntersectables.end(); ++ it)
	{
		DEL_PTR((*it).second);
	}

	// if not using kd tree root, delete tree (otherwise mKdTree has to do the job)
	if (!mCopyFromKdTree)
		DEL_PTR(mRoot);

	DEL_PTR(mSubdivisionCandidates);
	mSubdivisionStats.close();
}


void OspTree::ReadEnvironment()
{
	bool randomize = false;
	Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
	if (randomize)
		Randomize(); // initialise random generator for heuristics

	//-- termination criteria for autopartition
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxDepth", mTermMaxDepth);
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxLeaves", mTermMaxLeaves);
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.minObjects", mTermMinObjects);
	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minProbability", mTermMinProbability);
	
	Environment::GetSingleton()->GetIntValue("OspTree.Termination.missTolerance", mTermMissTolerance);
	
	//-- max cost ratio for early tree termination
	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.maxCostRatio", mTermMaxCostRatio);

	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minGlobalCostRatio",
		mTermMinGlobalCostRatio);
	

	//-- factors for bsp tree split plane heuristics
	Environment::GetSingleton()->GetFloatValue("OspTree.Termination.ct_div_ci", mCtDivCi);

	//-- partition criteria
	Environment::GetSingleton()->GetFloatValue("OspTree.Construction.epsilon", mEpsilon);

	// if only the driving axis is used for axis aligned split
	Environment::GetSingleton()->GetBoolValue("OspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
	Environment::GetSingleton()->GetFloatValue("OspTree.maxStaticMemory", mMaxMemory);
	Environment::GetSingleton()->GetBoolValue("OspTree.useCostHeuristics", mUseCostHeuristics);


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

	Environment::GetSingleton()->GetFloatValue("OspTree.Construction.splitBorder", mSplitBorder);
	Environment::GetSingleton()->GetFloatValue("OspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
	

	//-- debug output

	Debug << "******* OSP tree options ******** " << endl;

    Debug << "max depth: " << mTermMaxDepth << endl;
	Debug << "min probabiliy: " << mTermMinProbability << endl;
	Debug << "min objects: " << mTermMinObjects << endl;
	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	Debug << "miss tolerance: " << mTermMissTolerance << endl;
	Debug << "max leaves: " << mTermMaxLeaves << endl;

	Debug << "randomize: " << randomize << endl;

	Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
	Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
	Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
	Debug << "max memory: " << mMaxMemory << endl;
	Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
	Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
	
	Debug << "split borders: " << mSplitBorder << endl;
	Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
	Debug << endl;
}


void OspTree::SplitObjects(KdLeaf *parent,
						   const AxisAlignedPlane& splitPlane,
						   const ObjectContainer &objects,
						   ObjectContainer &front,
						   ObjectContainer &back)
{
	ObjectContainer::const_iterator oit, oit_end = objects.end();

    for (oit = objects.begin(); oit != oit_end; ++ oit) 
	{
		Intersectable *object = *oit;
		
		// initialise leaf references of objects
		if (parent->IsRoot())
		{
			object->mReferences = 1;
		}

		// remove parent ref
		-- object->mReferences;

		// determine the side of this ray with respect to the plane
		const AxisAlignedBox3 box = object->GetBox();

		if (box.Max(splitPlane.mAxis) >= splitPlane.mPosition)
		{
            front.push_back(object);
			++ object->mReferences;
		}

		if (box.Min(splitPlane.mAxis) < splitPlane.mPosition)
		{
			back.push_back(object);
			++ object->mReferences;
		}
	}

	mOspStats.objectRefs -= (int)objects.size();
	mOspStats.objectRefs += (int)(back.size() + front.size());
}


KdInterior *OspTree::SubdivideNode(
								   const AxisAlignedPlane &splitPlane,
								   const OspTraversalData &tData,
								   OspTraversalData &frontData,
								   OspTraversalData &backData
								   )
{
	KdLeaf *leaf = tData.mNode;
	
	// two new leaves
    mOspStats.nodes += 2;

	// add the new nodes to the tree
	KdInterior *node = new KdInterior(leaf->mParent);

	const int axis = splitPlane.mAxis;
	const float position = splitPlane.mPosition;

	node->mAxis = axis;
	node->mPosition = position;
	node->mBox = tData.mBoundingBox;


	//-- the front and back traversal data is filled with the new values

	// bounding boxes: split front and back node geometry
	tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition, 
		frontData.mBoundingBox, backData.mBoundingBox);

	frontData.mDepth = backData.mDepth = tData.mDepth + 1;
		
	//-- subdivide rays
	frontData.mRays = new RayInfoContainer();
	backData.mRays = new RayInfoContainer();

/*	SplitRays(splitPlane,
			  *tData.mRays,
			  *frontData.mRays,
			  *backData.mRays);
*/

	//-- classify objects

	int objectsBack = 0;
	int objectsFront = 0;

    ObjectContainer::const_iterator mi, mi_end = leaf->mObjects.end();

	for ( mi = leaf->mObjects.begin(); mi != mi_end; ++ mi) 
	{
		// determine the side of this ray with respect to the plane
		const AxisAlignedBox3 box = (*mi)->GetBox();
		
		if (box.Max(axis) >= position) 
			++ objectsFront;
    
		if (box.Min(axis) < position)
			++ objectsBack;
	}


	// TODO matt: compute pvs
	//frontData.mPvs = objectsFront;
	//backData.mPvs = objectsBack;

	//CheckViewCellsPvs(leaf, *tData.mRays);
	RemoveParentViewCellsPvs(leaf, *tData.mRays);


	/////////////
	//-- create front and back leaf

	KdLeaf *back = new KdLeaf(node, objectsBack);
	KdLeaf *front = new KdLeaf(node, objectsFront);

	KdInterior *parent = leaf->mParent;

	// replace a link from node's parent
	if (parent)
	{
		parent->ReplaceChildLink(leaf, node);
		node->mParent = parent;
	}
	else // new root
	{
		mRoot = node;
	}

	// and setup child links
	node->SetupChildLinks(back, front);

	SplitObjects(leaf, splitPlane, leaf->mObjects, front->mObjects, back->mObjects);
	
	FilterRays(front, *tData.mRays, *frontData.mRays);
	FilterRays(back, *tData.mRays, *backData.mRays);

	++ mOspStats.splits[splitPlane.mAxis];

	//-- eval view cells pvs
	
	// remove parent pvs update pvs of left and right leaves
	// Note: It is necessary to update first left and then right pvs.
	// We need to store the view cells seen by each object,
	// but also if the view cells are seen as part of two different
	// kd leaves, which is stored in the pdf component of the pvs entry. 
	// Because if an object is part of a view cell more than once,
	// it cannot be removed from the pvs by splitting a kd node where
	// the view cell sees only the other child of the node.
	// This is important during the subdivision

//MailablePvsData::NewMail();
	UpdateViewCellsPvs(front, *frontData.mRays);
	UpdateViewCellsPvs(back, *backData.mRays);

	////////////////////////////////////


	ProcessMultipleRefs(front);
    ProcessMultipleRefs(back);

	frontData.mNode = front;
	backData.mNode = back;


	// compute probability, i.e., volume of seen view cells
	frontData.mProbability = EvalViewCellsVolume(front, *frontData.mRays);
	backData.mProbability =  EvalViewCellsVolume(back, *backData.mRays);


	//delete leaf;
	return node;
}


KdNode *OspTree::Subdivide(SplitQueue &tQueue,
						   SubdivisionCandidate *splitCandidate,
						   const bool globalCriteriaMet)
{
	OspSubdivisionCandidate *sc = dynamic_cast<OspSubdivisionCandidate *>(splitCandidate);
	OspTraversalData &tData = sc->mParentData;

	KdNode *newNode = tData.mNode;

	if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
	{	
		OspTraversalData tFrontData;
		OspTraversalData tBackData;

		//-- continue subdivision
		
		// create new interior node and two leaf node
		const AxisAlignedPlane splitPlane = sc->mSplitPlane;
		
		newNode = SubdivideNode(splitPlane, 
								tData, 
								tFrontData, 
								tBackData);
	
		const int maxCostMisses = sc->mMaxCostMisses;

        // how often was max cost ratio missed in this branch?
		tFrontData.mMaxCostMisses = maxCostMisses;
		tBackData.mMaxCostMisses = maxCostMisses;
		
		mTotalCost -= sc->GetRenderCostDecrease();

		// subdivision statistics
		if (1) EvalSubdivisionStats(*sc);

		//-- push the new split candidates on the queue
		
		OspSubdivisionCandidate *frontCandidate = new OspSubdivisionCandidate(tFrontData);
		OspSubdivisionCandidate *backCandidate = new OspSubdivisionCandidate(tBackData);

		EvalSubdivisionCandidate(*frontCandidate);
		EvalSubdivisionCandidate(*backCandidate);

		tQueue.Push(frontCandidate);
		tQueue.Push(backCandidate);

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


	//-- terminate traversal
	if (newNode->IsLeaf())
	{
		EvaluateLeafStats(tData);
        
		const bool mStoreRays= true;

		//-- store additional info
		if (mStoreRays)
		{
			KdLeaf *leaf = dynamic_cast<KdLeaf *>(newNode);

			RayInfoContainer::const_iterator it, it_end = tData.mRays->end();

			for (it = tData.mRays->begin(); it != it_end; ++ it)
			{
				(*it).mRay->Ref();			
				leaf->mVssRays.push_back((*it).mRay);
			}
		}
	}
	
	//-- cleanup
	tData.Clear();
	
	return newNode;
}


void OspTree::EvalSubdivisionCandidate(OspSubdivisionCandidate &splitCandidate)
{
	// compute locally best split plane
	const float ratio = SelectSplitPlane(splitCandidate.mParentData, 
										 splitCandidate.mSplitPlane);

	const bool success = ratio < mTermMaxCostRatio;
	float oldRenderCost;

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

	splitCandidate.SetRenderCostDecrease(renderCostDecr);

#if 0
	const float priority = (float)-data.mDepth;
#else	
	// take render cost of node into account 
	// otherwise danger of being stuck in a local minimum!!
	const float factor = mRenderCostDecreaseWeight;
	const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;

#endif

	// compute global decrease in render cost
	splitCandidate.SetPriority(priority);

	splitCandidate.mMaxCostMisses = 
		success ? splitCandidate.mParentData.mMaxCostMisses : 
	splitCandidate.mParentData.mMaxCostMisses + 1;
}


inline bool OspTree::LocalTerminationCriteriaMet(const OspTraversalData &data) const
{
	// matt: TODO
	return ( 
		//(data.mNode->mObjects.size() < mTermMinObjects) ||
		//(data.mProbability <= mTermMinProbability) ||
		(data.mDepth >= mTermMaxDepth)
		 );
}


inline  bool OspTree::GlobalTerminationCriteriaMet(const OspTraversalData &data) const
{
	// matt: TODO
	return (
		(mOspStats.Leaves() >= mTermMaxLeaves)
		//mOutOfMemory || 
		//(mGlobalCostMisses >= mTermGlobalCostMissTolerance)
		);
}


void OspTree::EvaluateLeafStats(const OspTraversalData &data)
{
	// the node became a leaf -> evaluate stats for leafs
	KdLeaf *leaf = data.mNode;
	
	++ mCreatedLeaves;

	/*if (data.mPvs > mOspStats.maxPvs)
	{
		mOspStats.maxPvs = data.mPvs;
	}

	mOspStats.pvs += data.mPvs;

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

	if (data.mProbability <= mTermMinProbability)
		++ mOspStats.minProbabilityNodes;
	
	// accumulate depth to compute average depth
	mOspStats.accumDepth += data.mDepth;
 
	//	if ((int)(leaf->mObjects.size()) < mTermMinCost)
	//     ++ mOspStats.minCostNodes;
 
	if ((int)(leaf->mObjects.size()) > mOspStats.maxObjectRefs)
		mOspStats.maxObjectRefs = (int)leaf->mObjects.size();

#ifdef _DEBUG
	Debug << "BSP stats: "
		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
	  	  << "Prob: " << data.mProbability << " (min: " << mTermMinProbability << ")\n";
#endif
}


float OspTree::EvalLocalCostHeuristics(const OspTraversalData &tData,
									   const AxisAlignedBox3 &box,
									   const int axis,
									   float &position,
									   int &objectsFront,
									   int &objectsBack)
{
	RayInfoContainer usedRays;
	int mMaxTests = 500000; // HACK

	if (mMaxTests < (int)tData.mRays->size())
	{
		GetRayInfoSets(*tData.mRays, mMaxTests, usedRays);
	}
	else
	{
		usedRays = *tData.mRays;
	}

	// 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
	
	const float minBox = box.Min(axis);
	const float maxBox = box.Max(axis);
	Debug << "min: " << minBox << " max: " << maxBox << endl;

	const float sizeBox = maxBox - minBox;

	float minBand = minBox + mSplitBorder * (maxBox - minBox);
	float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);

	const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();

	//-- sort so we can use a sweep
	SortSubdivisionCandidates(tData, axis, minBand, maxBand);

	float totalVol = 0;
	float voll = 0;
	float volr = totalVol;

	ViewCellContainer touchedViewCells;
	
	const float totalRenderCost = PrepareHeuristics(tData, touchedViewCells);
	float renderCost = totalRenderCost;

	/// this is kind of a reverse pvs
	const int totalPvs = (int)tData.mNode->mObjects.size();
	
	int pvsl = 0;
	int pvsr = totalPvs;

	float sum = (float)totalVol * sizeBox;

	Debug << "here82 render cost: " << renderCost / viewSpaceVol << endl;

	/////////////////////////////////

	// note: initialised to take mid split 
	// if no good split can be found,
	position = minBox + 0.5f * sizeBox;
	
	// the relative cost ratio
	float ratio = 99999999.0f;
	bool splitPlaneFound = false;

	float volBack = voll;
	float volFront = volr;

	int pvsBack = pvsl;
	int pvsFront = pvsr;
	
	float minRenderCost = 1e20f;

	debugVol = 0;



	/////////////////////////////
	// the sweep heuristics

	
	//-- traverse through events and find best split plane
	
	vector<SortableEntry>::const_iterator ci, ci_end = mSubdivisionCandidates->end();
	//minRenderCost = RandomValue(0,viewSpaceVol);

	for (ci = mSubdivisionCandidates->begin(); ci != ci_end; ++ ci)
	{
		Intersectable *object = (*ci).mObject;

		EvalHeuristicsContribution(tData.mNode, 
								   *ci, 
								   renderCost,
								   touchedViewCells
								  );

		// Note: sufficient to compare size of bounding boxes of front and back side?
		if (((*ci).mPos >= minBand) && ((*ci).mPos <= maxBand))
		{
			float currentPos;
			Debug << "here29 : " << minRenderCost / viewSpaceVol << endl;

			// HACK: current positition is BETWEEN visibility events
			if (1 && /*((*ci).mType == SortableEntry::BOX_INTERSECT) &&*/ ((ci + 1) != ci_end))
				currentPos = ((*ci).mPos + (*(ci + 1)).mPos) * 0.5f;
			else
				currentPos = (*ci).mPos;			

			if (renderCost < minRenderCost)
			{
				splitPlaneFound = true;
				minRenderCost = renderCost;
				position = currentPos;
				Debug << "pos: " << position << endl;
			}
		}
	}
//splitPlaneFound = true;
	if (splitPlaneFound)
	{
		ratio = minRenderCost / totalRenderCost;
	}

	Debug << "\n§§§§ eval local cost §§§§" << endl
		  << "old rc: " << totalRenderCost / viewSpaceVol << " new rc: " << minRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << (totalRenderCost - minRenderCost) / viewSpaceVol  << endl;

	return ratio;
}


void OspTree::SortSubdivisionCandidates(const OspTraversalData &tData,
								  const int axis, 
								  float minBand, 
								  float maxBand)

{
	mSubdivisionCandidates->clear();

	RayInfoContainer *rays = tData.mRays;
	KdLeaf *leaf = tData.mNode;

	int requestedSize = 2 * (int)rays->size() + 2 * (int)leaf->mObjects.size();

	// creates a sorted split candidates array
	if (mSubdivisionCandidates->capacity() > 500000 &&
		requestedSize < (int)(mSubdivisionCandidates->capacity() / 10) )
	{
    	delete mSubdivisionCandidates;
   		mSubdivisionCandidates = new vector<SortableEntry>;
	}

	mSubdivisionCandidates->reserve(requestedSize);

	float pos;

	//-- insert ray queries
	//-- we are intersested only in rays which intersect an object that
	//-- is part of the kd node because they can induce a change in view cell
	//-- volume on the left and the right part. 
	//-- Note that they cannot induce a change in pvs size!!

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

	for (rit = rays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		const bool positive = ray->HasPosDir(axis);
	
		if (EndPointInsideNode(leaf, *ray, true))
		{
			pos = ray->mTermination[axis];

			mSubdivisionCandidates->push_back(
				SortableEntry(SortableEntry::BOX_INTERSECT, 
							  pos, 
							  ray->mTerminationObject, 
							  ray)
							  );
		}

		// if hitpoint with object is inside this node
		if (EndPointInsideNode(leaf, *ray, false))
		{
			pos = ray->mOrigin[axis];
	
			mSubdivisionCandidates->push_back(
				SortableEntry(SortableEntry::BOX_INTERSECT, 
							  pos, 
							  ray->mOriginObject, 
							  ray)
							  );
		}
	}


	//-- insert object queries
	//-- These queries can induce a change in pvs size.
	//-- Note that they cannot induce a change in view cell volume, as
	//-- there is no ray associated with these events!!

	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

    for (oit = leaf->mObjects.begin(); oit != leaf->mObjects.end(); ++ oit ) 
	{
		Intersectable *object = *oit;
		const AxisAlignedBox3 box = object->GetBox();

		mSubdivisionCandidates->push_back(
			SortableEntry(SortableEntry::BOX_MIN,
						  box.Min(axis),
						  object,
						  NULL)
						  );
   
    
		mSubdivisionCandidates->push_back(
			SortableEntry(SortableEntry::BOX_MAX,
						  box.Max(axis),
						  object,
						  NULL)
						  );
	}

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


const OspTreeStatistics &OspTree::GetStatistics() const
{
	return mOspStats;
}


void OspTree::PrepareHeuristics(const VssRay &ray, 
								ViewCellContainer &touchedViewCells)
{
	// collect view cells and set mail + counter
	ViewCellContainer viewCells;
	mVspTree->GetViewCells(ray, viewCells);
    	
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = (*vit);

		if (!vc->Mailed())
		{
			vc->Mail(); // set as belonging to front leaf
			vc->mCounter = 0;
			touchedViewCells.push_back(vc);
		}
		
		// view cell volume already added => just increase counter
		++ vc->mCounter;
	}
}


float OspTree::PrepareHeuristics(const OspTraversalData &tData,
								 ViewCellContainer &touchedViewCells)
{	
	float renderCost = 0;

	Intersectable::NewMail(3);
	ViewCell::NewMail(3);
	MailablePvsData::NewMail();

	KdLeaf *leaf = tData.mNode;

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();

	for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
     	PrepareHeuristics(*ray, touchedViewCells);
	}

	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		obj->Mail(); // set as belonging to front leaf
		
		ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

			if (!vdata) // should never happen
				continue;

			if (!vdata->Mailed())
			{
				vdata->Mail();
				renderCost += vc->GetVolume();
			}
		}
	}

	Debug << "here32 " << touchedViewCells.size() << endl;
	return renderCost;
}


void OspTree::AddObjectContribution(KdLeaf *leaf,
						  Intersectable *obj,
						  ViewCellContainer &touchedViewCells,
						  float &renderCost)
{
	//Debug << "add contri to obj: " << obj->mMailbox - Intersectable::sMailId << endl;
	obj->Mail(2); // set as belonging to both leafs

	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();
	
	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;
	
		// if obj not previously associated with this view cell => increase render cost
		if (vc->Mailed(1) && !ViewCellHasMultipleReferences(obj, vc, false))
		{	
			//Debug << "add to rendercost: " << vc->GetVolume() / mVspTree->GetBoundingBox().GetVolume() << endl;
			renderCost += vc->GetVolume();
		}
	}
}


void OspTree::SubtractObjectContribution(KdLeaf *leaf,
								Intersectable * obj,
								ViewCellContainer &touchedViewCells,
								float &renderCost)
{
	//Debug << "subtract contri from obj: " << obj->mMailbox - Intersectable::sMailId << endl;
    obj->Mail(1); // set as belonging to back leaf
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();
    	
	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;

		// if obj was previously associated with this view cell but is not now 
		// => decrease render cost
		if (vc->Mailed() &&	!ViewCellHasMultipleReferences(obj, vc, false))
		{
			//Debug << "subtract from rendercost: " << vc->GetVolume() / mVspTree->GetBoundingBox().GetVolume() << endl;
			renderCost -= vc->GetVolume();
		}
	}
}


void OspTree::EvalRayContribution(KdLeaf *leaf,
								  const VssRay &ray, 
								  float &renderCost)
{
	ViewCellContainer viewCells;

	mVspTree->GetViewCells(ray, viewCells);

	// classify view cells and compute volume contri accordingly
	// possible view cell classifications:
	// view cell mailed => view cell can be seen from left child node
	// view cell counter > 0 view cell can be seen from right child node
	// combined: view cell volume belongs to both nodes
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
	
	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
		EvalViewCellContribution(leaf, *vit, renderCost);
	}
}


void OspTree::EvalViewCellContribution(KdLeaf *leaf, 
									   ViewCell *viewCell,
									   float &renderCost)
{
	//Debug << "**************" << endl;
	//Debug << "vc contri: " << renderCost / mVspTree->GetBoundingBox().GetVolume() << " " << viewCell->mMailbox - ViewCell::sMailId << " counter: " << viewCell->mCounter << endl;

	const float vol = viewCell->GetVolume();

	if (viewCell->Mailed())
	{
		viewCell->Mail(2); // view cell can be seen from both nodes

        // we now see view cell from both nodes => add contri
		ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

		for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
		{
			Intersectable *obj = *oit;

			// was render cost already added?
			if (!ViewCellHasMultipleReferences(obj, viewCell, false) &&
				obj->Mailed(1))
			{
				//Debug << "obj mail 1!! " << vol / mVspTree->GetBoundingBox().GetVolume() << endl;
				renderCost += vol;
			}
		}
	}

	if (-- viewCell->mCounter == 0)
	{
		viewCell->Mail(1); // view cell can be seen from back node only

		//MailablePvsData::NewMail();
		ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

		for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
		{
			Intersectable *obj = *oit;

			// can render cost be be reduced?
			if (!ViewCellHasMultipleReferences(obj, viewCell, false) &&	obj->Mailed())
			{
				renderCost -= vol;
			}
		}
	}
	//Debug << "new rc: " << renderCost / mVspTree->GetBoundingBox().GetVolume() << " counter2: " << viewCell->mCounter << endl;
	//Debug << "**************"<<endl;
}


void OspTree::EvalHeuristicsContribution(KdLeaf *leaf,
										 const SortableEntry &ci,
										 float &renderCost,
										 ViewCellContainer &touchedViewCells)
{
	Intersectable *obj = ci.mObject;
	VssRay *ray = ci.mRay;

	// switch between different types of events
	switch (ci.mType) 
	{
		case SortableEntry::BOX_MIN:
			cout << "<";
			AddObjectContribution(leaf, ci.mObject, touchedViewCells, renderCost);	
			break;
			
		case SortableEntry::BOX_MAX:
			cout << ">";
			SubtractObjectContribution(leaf, ci.mObject, touchedViewCells, renderCost);
			break;

		// compute volume contribution from view cells
		case SortableEntry::BOX_INTERSECT:
			cout << "i";
			// process ray if the hit point with termination / origin object 
			// is inside this kd leaf
			EvalRayContribution(leaf, *ray, renderCost);
			break;

		default:
			cout << "should not come here" << endl;
			break;
	}

	//cout << "vol left: " << volLeft << " vol right " << volRight << endl;
}


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


AxisAlignedBox3 OspTree::GetBoundingBox() const
{
	return mBoundingBox;
}


float OspTree::SelectSplitPlane(const OspTraversalData &tData,
								AxisAlignedPlane &plane)
{
	float nPosition[3];
	float nCostRatio[3];

	// create bounding box of node geometry
	AxisAlignedBox3 box = tData.mBoundingBox;
		
	int sAxis = 0;
	int bestAxis = -1;

	if (mOnlyDrivingAxis)
	{
		sAxis = box.Size().DrivingAxis();
	}
	

	// -- evaluate split cost for all three axis
	for (int axis = 0; axis < 3; ++ axis)
	{
		if (!mOnlyDrivingAxis || (axis == sAxis))
		{
			if (1 || mUseCostHeuristics)
			{
				//-- place split plane using heuristics
				int objectsFront, objectsBack;

				nCostRatio[axis] =
					EvalLocalCostHeuristics(tData,
									   tData.mBoundingBox,
									   axis,
									   nPosition[axis],
									   objectsFront,
									   objectsBack);			
			}
			/*	else 
			{
				//-- split plane position is spatial median

				nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;

				nCostRatio[axis] = EvalLocalSplitCost(tData,
													  box,
													  axis,
													  nPosition[axis],
													  nProbFront[axis], 
													  nProbBack[axis]);
			}*/
						
			if (bestAxis == -1)
			{
				bestAxis = axis;
			}
			else if (nCostRatio[axis] < nCostRatio[bestAxis])
			{
				bestAxis = axis;
			} 
		}
	}

	//-- assign values

	plane.mAxis = bestAxis;
	// split plane position
    plane.mPosition = nPosition[bestAxis];

	Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
	return nCostRatio[bestAxis];
}


bool OspTree::EndPointInsideNode(KdLeaf *leaf, 
								 VssRay &ray, 
								 bool isTermination) const
{
	// test if the leaf where the hitpoint is located is the current leaf
	if (isTermination)
	{
		 return ray.mTerminationObject && (GetLeaf(ray.mTermination, ray.mTerminationNode) == leaf);
	}
	else
	{
		return false; // no origin object!
		return ray.mOriginObject && (GetLeaf(ray.mOrigin, ray.mOriginNode) == leaf);
	}
}


void OspTree::EvalSubdivisionStats(const SubdivisionCandidate &sc)
{
	const float costDecr = sc.GetRenderCostDecrease();

	AddSubdivisionStats(mOspStats.Leaves(),
						costDecr,
						mTotalCost
						);
}


void OspTree::AddSubdivisionStats(const int leaves,
								  const float renderCostDecr,
								  const float totalRenderCost)
{
	mSubdivisionStats 
			<< "#Leaves\n" << leaves << endl
			<< "#RenderCostDecrease\n" << renderCostDecr << endl 
			<< "#TotalRenderCost\n" << totalRenderCost << endl;
			//<< "#AvgRenderCost\n" << avgRenderCost << endl;
}


int OspTree::ClassifyRay(VssRay *ray, 
						 KdLeaf *leaf, 
						 const AxisAlignedPlane &plane) const
{
	const bool originInside = EndPointInsideNode(leaf, *ray, false);
    const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

	const bool originGt = (ray->mOrigin[plane.mAxis] >= plane.mPosition);
	const bool originLt = (ray->mOrigin[plane.mAxis] <= plane.mPosition);

	const bool terminationGt = (ray->mTermination[plane.mAxis] >= plane.mPosition);
	const bool terminationLt = (ray->mTermination[plane.mAxis] <= plane.mPosition);

	// add view cell volume to front volume
	const bool addToFront = 
		((originInside && originGt) || (terminationInside && terminationGt));

	// add view cell volume to back volume
	const bool addToBack = 
		((originInside && originLt/*!originGt*/) || (terminationInside && terminationLt/*terminationGt*/));

	// classify ray with respect to the child nodes the view cells contribute
	if (addToFront && addToBack)
		return 0;
	else if (addToBack)
		return -1;
	else if (addToFront)
		return 1;

	return -2;
}


int OspTree::ClassifyRays(const RayInfoContainer &rays, 
						  KdLeaf *leaf, 
						  const AxisAlignedPlane &plane,
						  RayInfoContainer &frontRays,
						  RayInfoContainer &backRays) const
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		bool originGt = false;
		bool terminationGt = false;
		
		Intersectable *obj = ray->mOriginObject;

		if (0 && obj)
		{
			originGt = ray->mOrigin[plane.mAxis] >= plane.mPosition;
		}

		obj = ray->mTerminationObject;
		
		if (obj)
		{
			terminationGt = ray->mTermination[plane.mAxis] >= plane.mPosition;
		}

		// either begin or end point on front side
		const bool addToFront = originGt || terminationGt;

		// add view cell volume to back volume
		const bool addToBack = !originGt || !terminationGt;
	
		// classify ray with respect to the child nodes the view cells contribute
		if (addToFront)
			frontRays.push_back(*rit);
		if (addToBack)
			backRays.push_back(*rit);
	}

	return 0;
}


#if 0

float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
									  const OspTraversalData &tData,
									  float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;

	// probability that view point lies in back / front node
	float pOverall = 0;//data.mProbability; // todo matt§: value ok?
	float pFront = 0;
	float pBack = 0;
	float pFrontAndBack = 0;

	const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();

	Intersectable::NewMail(3);
	
	KdLeaf *leaf = tData.mNode;
	const int totalPvs = (int)leaf->mObjects.size();
	
	RayInfoContainer frontRays, backRays;
	ClassifyRays(*tData.mRays, leaf, candidatePlane, frontRays, backRays);

	ViewCellContainer touchedViewCells;
	
	// sum up volume seen from the objects of left and right children
	// => the volume is the weight for the render cost equation
	ViewCell::NewMail(3);

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();

	for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// add volume to volumes of left and / or right children
		// if one of the ray end points is inside
		const int classification = ClassifyRay(ray, leaf, candidatePlane);

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);
			
		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// if not previously mailed
			if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
			{
				touchedViewCells.push_back(vc);
			}

			// classify / mail the view cell
			MailViewCell(vc, classification);	
		}
	}

	// evaluate view cells volume contribution 
	// with respect to the mail box classification
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
	}

	//////////////////////////////////////////////
	//
	// evaluate contribution of objects which are part of other kd nodes
	// which are seen by the same view cell
	// These contributions cannot be reduced by splitting, because
	// the object would still be part of the pvs

	float additionalFrontRenderCost = 0;
	float additionalBackRenderCost = 0;
	
	MailablePvsData::NewMail();

	for (rit = tData.mRays->begin(); rit != tData.mRays->end(); ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// traverse through view cells
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *viewCell = *vit;

			// for a view cell:
			// if object is also in the kd pvs entries from other kd leaves,
			// render cost cannot be reduced for this view cell
			// => the render cost was falsly reduced, add them again

			if (ray->mTerminationObject)
			{	
				Intersectable *obj = ray->mTerminationObject;
				
				const bool renderCostWrong = 
					ViewCellHasMultipleReferences(obj, viewCell, true);
			
				// if there is already an entry of this object in the view cell pvs 
				if (renderCostWrong)
				{
					// view cell was counted only for front or back => correct tjos
					if (viewCell->Mailed(1) && obj->Mailed())
					{
						additionalFrontRenderCost += viewCell->GetVolume();
					}
					else if (viewCell->Mailed() && obj->Mailed(1))
					{
						additionalBackRenderCost += viewCell->GetVolume();
					}
				}
			}
	
			if (ray->mOriginObject)
			{
				Intersectable *obj = ray->mOriginObject;

				const bool renderCostWrong = 
					ViewCellHasMultipleReferences(obj, viewCell, true);
				// if there is already an entry of this object in the view cell pvs 
				if (renderCostWrong)
				{
					if (viewCell->Mailed(1) && obj->Mailed())
					{
						additionalFrontRenderCost += viewCell->GetVolume();
					}
					else if (viewCell->Mailed() && obj->Mailed(1))
					{
						additionalBackRenderCost += viewCell->GetVolume();
					}
				}
			}
		}
	}

	/////////////////////////////


	// these are non-overlapping sets
	pOverall = pFront + pBack + pFrontAndBack;


	//-- pvs rendering heuristics

	const float oldRenderCost = pOverall * totalPvs;
	
	// sum up the terms: 
	// the view cells seeing
	// a) the left node are weighted by the #left node objects
	// b) the right node are weighted by the #right node objects
	// c) both nodes are weighted by the #parent node objects
	const float newRenderCost = pvsFront * pFront + pvsBack * pBack + totalPvs * pFrontAndBack
		+ additionalFrontRenderCost + additionalBackRenderCost;

	// normalize volume with view space volume
	const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;

	Debug << "\n(((( eval render cost decrease ))))" << endl
		  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl 
		  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol 
		  << " front and back p " << pFrontAndBack / viewSpaceVol << " p: " << tData.mProbability / viewSpaceVol << endl
		  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << renderCostDecrease << endl
		  << "additional front " << additionalFrontRenderCost / viewSpaceVol 
		  << " additional back " << additionalBackRenderCost  / viewSpaceVol << endl;

	if (oldRenderCost < newRenderCost * 0.99)
		Debug <<"error!!"<<endl;

	//if ((((pOverall - tData.mProbability) / viewSpaceVol) > 0.00001))
	//	Debug << "ERROR!!"<<endl;

	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

		
	return renderCostDecrease;
}

#else

float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
									  const OspTraversalData &tData,
									  float &normalizedOldRenderCost) const
{
	float pvsFront = 0;
	float pvsBack = 0;

	// probability that view point lies in back / front node
	float pOverall = 0;//data.mProbability; // todo matt§: value ok?
	float pFront = 0;
	float pBack = 0;
	float pFrontAndBack = 0;

	Intersectable::NewMail(3);
	KdLeaf *leaf = tData.mNode;

	const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();
	const int totalPvs = (int)leaf->mObjects.size();
	
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	// evaluate reverse pvs and view cell volume on left and right cell
	// note: should I take all leaf objects or rather the objects hit by rays?
	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		int pvsContri = 0;

		Intersectable *obj = *oit;
		const AxisAlignedBox3 box = obj->GetBox();

		// test if box falls in left / right child node
		if (box.Max(candidatePlane.mAxis) >= candidatePlane.mPosition)
		{
			++ pvsFront;
			obj->Mail();		
		}

		if (box.Min(candidatePlane.mAxis) < candidatePlane.mPosition)
		{
			++ pvsBack;

			if (obj->Mailed())
				obj->Mail(2);
			else 
				obj->Mail(1);
		}
	}

	
	ViewCellContainer touchedViewCells;

	// sum up volume seen from the objects of left and right children
	// => the volume is the weight for the render cost equation
	ViewCell::NewMail(3);

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
	//map<Intersectable *, ViewCellContainer> objectsMap;

	for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// add volume to volumes of left and / or right children
		// if one of the ray end points is inside
		const int classification = ClassifyRay(ray, leaf, candidatePlane);

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);
			
		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// if not previously mailed
			if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
			{
				touchedViewCells.push_back(vc);
			}

			// classify / mail the view cell
			MailViewCell(vc, classification);
		}
	}

	// evaluate view cells volume contribution 
	// with respect to the mail box classification
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
	}

	//////////////////////////////////////////////
	//
	// evaluate contribution of objects which are part of other kd nodes
	// which are seen by the same view cell
	// These contributions cannot be reduced by splitting, because
	// the object would still be part of the pvs

	float newRc = 0;
	float rc = 0;

	MailablePvsData::NewMail();
	//ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			
			MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

			if (!vdata)
			{
				Debug << "here86 error!!"<<endl;
				continue;
			}

			if (!vdata->Mailed())
			{
				vdata->Mail();
				rc += vc->GetVolume();

				if ((vdata->mSumPdf > 1.5) || 
					vc->Mailed(2) || 
					obj->Mailed(2) ||
					(obj->Mailed() && vc->Mailed()) ||
					(obj->Mailed(1) && vc->Mailed(1))
					)
				{
					newRc += vc->GetVolume();
				}
			}
		}
	}


	/////////////////////////////


	// these are non-overlapping sets
	pOverall = pFront + pBack + pFrontAndBack;


	//-- pvs rendering heuristics

	const float oldRenderCost = rc;//pOverall * totalPvs;
	
	// sum up the terms: 
	// the view cells seeing
	// a) the left node are weighted by the #left node objects
	// b) the right node are weighted by the #right node objects
	// c) both nodes are weighted by the #parent node objects
	const float newRenderCost = newRc//pvsFront * pFront + pvsBack * pBack + totalPvs * pFrontAndBack
;//		+ additionalFrontRenderCost + additionalBackRenderCost;

	// normalize volume with view space volume
	const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;

	Debug << "\n(((( eval render cost decrease ))))" << endl
		  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl 
		  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol 
		  << " front and back p " << pFrontAndBack / viewSpaceVol << " p: " << tData.mProbability / viewSpaceVol << endl
		  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
		  << "render cost decrease: " << renderCostDecrease << endl;
		  //<< "additional front " << additionalFrontRenderCost / viewSpaceVol 
		  //<< " additional back " << additionalBackRenderCost  / viewSpaceVol << endl;

	if (oldRenderCost < newRenderCost * 0.99)
		Debug <<"error!!"<<endl;

	//if ((((pOverall - tData.mProbability) / viewSpaceVol) > 0.00001))
	//	Debug << "ERROR!!"<<endl;

	normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

		
	return renderCostDecrease;
}
#endif


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

		ObjectContainer::const_iterator oit, oit_end = objects.end();

		//-- compute bounding box
        for (oit = objects.begin(); oit != oit_end; ++ oit)
		{
			Intersectable *obj = *oit;

			// compute bounding box of view space
			mBoundingBox.Include(obj->GetBox());
		}
	}
}


void OspTree::ProcessMultipleRefs(KdLeaf *leaf) const
{
	// find objects from multiple kd-leaves
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *object = *oit;
		
		if (object->mReferences > 1)
		{
			leaf->mMultipleObjects.push_back(object);
		}
	}
}


void OspTree::CollectLeaves(vector<KdLeaf *> &leaves) const
{
	stack<KdNode *> nodeStack;
	nodeStack.push(mRoot);

	while (!nodeStack.empty()) 
	{
		KdNode *node = nodeStack.top();
		nodeStack.pop();
		if (node->IsLeaf()) 
		{
			KdLeaf *leaf = (KdLeaf *)node;
			leaves.push_back(leaf);
		}
		else 
		{
			KdInterior *interior = (KdInterior *)node;
			nodeStack.push(interior->mBack);
			nodeStack.push(interior->mFront);
		}
	}
}


AxisAlignedBox3 OspTree::GetBoundingBox(KdNode *node) const
{
	if (!node->mParent)
		return mBoundingBox;

	if (!node->IsLeaf())
	{
		return (dynamic_cast<KdInterior *>(node))->mBox;		
	}

	KdInterior *parent = dynamic_cast<KdInterior *>(node->mParent);

	AxisAlignedBox3 box(parent->mBox);

	if (parent->mFront == node)
		box.SetMin(parent->mAxis, parent->mPosition);
    else
		box.SetMax(parent->mAxis, parent->mPosition);

	return box;
}


void OspTree::CollectViewCells(KdLeaf *leaf, ViewCellContainer &viewCells)
{
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	ViewCell::NewMail();

	// loop through all object and collect view cell pvs of this node
	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

		ViewCellPvsMap::const_iterator vit, vit_end = obj->mViewCellPvs.mEntries.end();

		for (vit = obj->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
		{
            ViewCell *vc = (*vit).first;

			if (!vc->Mailed())
			{
				vc->Mail();
				viewCells.push_back(vc);
			}
		}
	}
}


void OspTree::CollectDirtyCandidates(OspSubdivisionCandidate *sc, 
									 vector<SubdivisionCandidate *> &dirtyList)
{
	OspTraversalData &tData = sc->mParentData;
	KdLeaf *node = tData.mNode;
	
	ViewCell::NewMail();
	ViewCellContainer viewCells;
	ViewCellContainer tmpViewCells;

	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();

	// find all view cells associated with the rays, add them to view cells
	for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
		mVspTree->GetViewCells(*ray, tmpViewCells);

		ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();

		for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
		{
			VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);

			if (!vc->Mailed())
			{
				vc->Mail();
				viewCells.push_back(vc);
			}
		}
	}


	// split candidates holding this view cells 
	// are thrown into dirty list
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
		VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);

		VspLeaf *leaf = vc->mLeaf;
		dirtyList.push_back(leaf->GetSubdivisionCandidate());
	}
}


KdNode *OspTree::GetRoot() const
{
	return mRoot;
}


KdLeaf *OspTree::GetLeaf(const Vector3 &pt, KdNode *node) const
{
	// start from root of tree
	if (node == NULL)
	{
		node = mRoot;
	}

	stack<KdNode *> nodeStack;
	nodeStack.push(node);
  
	KdLeaf *leaf = NULL;
  
	while (!nodeStack.empty())  
	{
		KdNode *node = nodeStack.top();
		nodeStack.pop();
	
		if (node->IsLeaf()) 
		{
			leaf = dynamic_cast<KdLeaf *>(node);
		} 
		else 	
		{	
			// find point
			KdInterior *interior = dynamic_cast<KdInterior *>(node);
        
			if (interior->mPosition > pt[interior->mAxis])
			{
				nodeStack.push(interior->mBack);
			}
			else
			{
				nodeStack.push(interior->mFront);
			}
		}
	}
  
	return leaf;
}


void OspTree::PreprocessRays(KdLeaf *root,
							 const VssRayContainer &sampleRays,
							 RayInfoContainer &rays)
{
	VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
	RayInfoContainer clippedRays;

	long startTime = GetTime();

	cout << "storing rays ... ";

	Intersectable::NewMail();

	//-- store rays and objects
	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = *rit;

		float minT, maxT;
		static Ray hray;

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

			clippedRays.push_back(RayInfo(ray, minT / len, maxT / len));

			// HACK: reset nodes for the case we have used object kd tree for
			// tree building.
			ray->mOriginNode = NULL;
			ray->mTerminationNode = NULL;
		}
	}

	FilterRays(root, clippedRays, rays);

	cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
}



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

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

	for (rit = rays.begin(); rit != rit_end; ++ rit)
	{
		RayInfo bRay = *rit;
		
		VssRay *ray = bRay.mRay;
		float t;

		// get classification and receive new t
		//-- test if start point behind or in front of plane
		const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
			
		if (side == 0)
		{
			++ splits;

			if (ray->HasPosDir(plane.mAxis))
			{
				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()));
			}
		}
		else if (side == 1)
		{
			frontRays.push_back(bRay);
		}
		else
		{
			backRays.push_back(bRay);
		}
	}

	return splits;
}


int OspTree::FilterRays(KdLeaf *leaf, 
						const RayInfoContainer &rays,
						RayInfoContainer &filteredRays)
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	for (rit = rays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// test if intersection point with one of the objects is inside this node
		const bool originInside = EndPointInsideNode(leaf, *ray, false);
        const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

		if (originInside || terminationInside)
		{
			filteredRays.push_back(ray);
		}
	}

	return 0;
}
					 

int OspTree::CheckViewCellsPvs(KdLeaf *leaf, 
							   const RayInfoContainer &rays) const
{
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	Intersectable::NewMail();
	ObjectContainer touchedObjects;
	

	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		if (ray->mTerminationObject)
		{
			Intersectable *obj = ray->mTerminationObject;

			if (!obj->Mailed())
			{
				obj->Mail();
				touchedObjects.push_back(obj);
			}
		}

		if (ray->mOriginObject)
		{
			Intersectable *obj = ray->mOriginObject;

			if (!obj->Mailed())
			{
				obj->Mail();
				touchedObjects.push_back(obj);
			}
		}
	}
	Debug << "here65 " << touchedObjects.size() << endl;
	ObjectContainer::const_iterator it, it_end = touchedObjects.end();
	for (it = touchedObjects.begin(); it != it_end; ++ it)
	{
		Debug << "\nhere94 obj: " << (*it) << " size: " << (*it)->mViewCellPvs.GetSize() << endl << endl;
		ViewCellPvsMap::const_iterator mit, mit_end = (*it)->mViewCellPvs.mEntries.end();

		for (mit = (*it)->mViewCellPvs.mEntries.begin(); mit != mit_end; ++ mit)
		{
			Debug << "newsumpdf: " << (*mit).second.mSumPdf << endl;
		}
	}

	return 0;
}


KdIntersectable *OspTree::GetOrCreateKdIntersectable(KdNode *node)
{
	// search nodes
	std::map<KdNode *, KdIntersectable *>::
		const_iterator it = mKdIntersectables.find(node);

	if (it != mKdIntersectables.end()) 
	{
		return (*it).second;
	}

	// not in map => create new entry
	KdIntersectable *kdObj = new KdIntersectable(node);
	mKdIntersectables[node] = kdObj;

	return kdObj;
}


/*
void OspTree::AddViewCellVolume(ViewCell *vc,
								const int cf,
								float &frontVol,
								float &backVol,
								float &totalVol) const
{
	//const float renderCost = mViewCellsManager->SimpleRay &raynderCost(obj);
	const float vol = vc->GetVolume();

	// view cell not found yet => new
	if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
	{
		totalVol += vol;
	}

	if (cf >= 0) // front volume
	{
		if (!vc->Mailed() && !vc->Mailed(2))
		{
			frontVol += vol;
		
			// already in back volume => in both volumes
			if (vc->Mailed(1))
				vc->Mail(2);
			else
				vc->Mail();
		}
	}

	if (cf <= 0) // back volume
	{
		if (!vc->Mailed(1) && !vc->Mailed(2))
		{
			backVol += vol;
		
			// already in front volume => in both volume
			if (vc->Mailed())
				vc->Mail(2);
			else
				vc->Mail(1);
		}
	}
}
*/


void OspTree::MailViewCell(ViewCell *vc, const int cf) const
{
	if (vc->Mailed(2)) // already classified
		return;

	if (cf >= 0) // front volume
	{
		// already in back volume => in both volumes
		if (vc->Mailed(1))
		{
			vc->Mail(2);
		}
		else
		{
			vc->Mail();
		}
	}

	if (cf <= 0) // back volume
	{
		// already in front volume => in both volume
		if (vc->Mailed())
		{
			vc->Mail(2);
		}
		else
		{
			vc->Mail(1);
		}
	}
}


void OspTree::AddViewCellVolumeContri(ViewCell *vc,
									  float &frontVol,
									  float &backVol,
									  float &frontAndBackVol) const
{
	if (vc->Mailed())
	{
		frontVol += vc->GetVolume();
	}
	else if (vc->Mailed(1))
	{
		backVol += vc->GetVolume();
	}
	else if (vc->Mailed(2))
	{
		frontAndBackVol += vc->GetVolume();
	}
}


float OspTree::EvalViewCellsVolume(KdLeaf *leaf, 
								   const RayInfoContainer &rays) const
{
	float vol = 0;
	ViewCell::NewMail();

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

	for (rit = rays.begin(); rit != rays.end(); ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
			
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			if (!vc->Mailed())
			{
				vc->Mail();
				vol += vc->GetVolume();
			}
		}		
	}

	return vol;
}


bool OspTree::AddViewCellToObjectPvs(Intersectable *obj, 
									 ViewCell *vc, 
									 float &contribution,
									 bool onlyMailed) const
{	
	contribution = 0; // todo

	MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

	if (!vdata)
	{
		vdata = obj->mViewCellPvs.AddSample2(vc, 1);
	}
	else if (!onlyMailed || !vdata->Mailed())
	{
		obj->mViewCellPvs.AddSample(vc, 1);
	}
	
	vdata->Mail();

	return true;
}


void OspTree::CollectTouchedViewCells(const RayInfoContainer &rays, 
									  ViewCellContainer &touchedViewCells) const
{
	ViewCell::NewMail();
	
	RayInfoContainer::const_iterator rit, rit_end = rays.end();

	for (rit = rays.begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			if (!vc->Mailed())
			{
				vc->Mail();
				touchedViewCells.push_back(vc);
			}
		}
	}
}


int OspTree::UpdateViewCellsPvs(KdLeaf *leaf, 
								const RayInfoContainer &rays) const

{
	MailablePvsData::NewMail();
	
	/*RayInfoContainer::const_iterator rit, rit_end = rays.end();

	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			Intersectable *obj = ray->mTerminationObject;
			if (obj)
			{
				float contri;
				AddViewCellToObjectPvs(obj, vc, contri, true);
			}

			obj = ray->mOriginObject;
			if (obj)
			{
				float contri;
				AddViewCellToObjectPvs(obj, vc, contri, true);
			}
		}
	}*/
	
	ViewCellContainer touchedViewCells;
	CollectTouchedViewCells(rays, touchedViewCells);

	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
	{
		Intersectable *obj = *oit;
		ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			float contri;
			AddViewCellToObjectPvs(obj, vc, contri, true);
		}
	}

	return 0;
}


int OspTree::RemoveParentViewCellsPvs(KdLeaf *leaf, 
									  const RayInfoContainer &rays
									  ) const

{
	MailablePvsData::NewMail();

	/*RayInfoContainer::const_iterator rit, rit_end = rays.end();
	for (rit = rays.begin(); rit < rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;

		// test if intersection point with one of the objects is inside this node
		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);

		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			Intersectable *obj = ray->mTerminationObject;

			if (obj)
			{
				MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

				if (vdata && !vdata->Mailed())
				{
					vdata->Mail();
					obj->mViewCellPvs.RemoveSample(vc, 1);
				}
			}

			obj = ray->mOriginObject;

			if (obj)
			{
				MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

				if (vdata && !vdata->Mailed())
				{
					vdata->Mail();
					obj->mViewCellPvs.RemoveSample(vc, 1);
				}
			}
		}
	}*/

	ViewCellContainer touchedViewCells;
	CollectTouchedViewCells(rays, touchedViewCells);

	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

		// traverse through view cells and classify them according
		// to them being seen from to back / front / front and back node
        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			
			MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

			if (vdata && !vdata->Mailed())
			{
				vdata->Mail();
				obj->mViewCellPvs.RemoveSample(vc, 1);
			}
		}
	}

	return 0;
}


float OspTree::EvalRenderCost(const VssRayContainer &myrays)
{
	float rc = 0;
	float prop = mVspTree->GetBoundingBox().GetVolume();	

	KdLeafContainer leaves;
	CollectLeaves(leaves);

	ViewCellContainer vcleaves;
	mVspTree->CollectViewCells(vcleaves, false);
	

	ViewCellContainer::const_iterator vit, vit_end = vcleaves.end();

	for (vit = vcleaves.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;
		vc->GetPvs().Clear();
	}

	KdLeafContainer::const_iterator kit, kit_end = leaves.end();

	MailablePvsData::NewMail();

	for (kit = leaves.begin(); kit != kit_end; ++ kit)
	{
		KdLeaf *leaf = *kit;
		float newCost = 0;
		float vol = 0;

		ViewCell::NewMail();
		VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
		ViewCellContainer touchedViewCells;

		for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
		{
			VssRay *ray = *rit;
			
			// test if intersection point with one of the objects is inside this node
			const bool originInside = EndPointInsideNode(leaf, *ray, false);
			const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

			if (originInside || terminationInside)
			{
				ViewCellContainer viewCells;
				mVspTree->GetViewCells(*ray, viewCells);
			
				ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

				for (vit = viewCells.begin(); vit != vit_end; ++ vit)
				{
					ViewCell *vc = *vit;

					// if not previously mailed
					if (!vc->Mailed())
					{
						vc->Mail();
						touchedViewCells.push_back(vc);
					}

        		}
			}
		}
		
		ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

		for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
		{
			Intersectable *obj = *oit;
			ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

			for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
			{
				ViewCell *vc = *vit;
				PvsData *vdata = vc->GetPvs().Find(obj);

				if (!vdata)
				{
					vc->GetPvs().AddSample(obj, 1);
					newCost += vc->GetVolume();
				}

/*				MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

				if (!vdata || !vdata->Mailed())
				{
					newCost += vc->GetVolume();
					if (!vdata)
						vdata = obj->mViewCellPvs.AddSample2(vc, 1);
					vdata->Mail();
				}*/
			}
		}

		rc += newCost;
	}

	return rc / prop;
}


float OspTree::EvalLeafCost(const OspTraversalData &tData)
{
	KdLeaf *leaf = tData.mNode;

	float rc = 0;
	//float prop = mVspTree->GetBoundingBox().GetVolume();	
	//float vol = 0;
	//ViewCell::NewMail();
	
	RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
	ViewCellContainer touchedViewCells;

	for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
	{
		VssRay *ray = (*rit).mRay;
			
		// test if intersection point with one of the objects is inside this node

		ViewCellContainer viewCells;
		mVspTree->GetViewCells(*ray, viewCells);
			
		ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

		for (vit = viewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;

			// if not previously mailed
			if (!vc->Mailed())
			{
				vc->Mail();
				touchedViewCells.push_back(vc);
			}
		}
	}

	// clear pvs of involved view cells
	ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

	for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
	{
		ViewCell *vc = *vit;
		vc->GetPvs().Clear();
	}
	
	ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

	for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
	{
		Intersectable *obj = *oit;

		for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
		{
			ViewCell *vc = *vit;
			PvsData *vdata = vc->GetPvs().Find(obj);

			if (!vdata)
			{
				vc->GetPvs().AddSample(obj, 1);
				rc += vc->GetVolume();
				//prop += vc->GetVolume();
			}
		}
	}
	
	return rc;
}


bool OspTree::Export(OUT_STREAM &stream)
{
	ExportNode(mRoot, stream);

	return true;
}


void OspTree::ExportNode(KdNode *node, OUT_STREAM &stream)
{
	if (node->IsLeaf())
	{
		KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);

		stream << "<Leaf ";
		stream << "objects=\"";
		
		// export objects in kd leaves
		//if (mExportObjects) ExportObjects(leaf, stream);
		
		stream << "\" />" << endl;
	}
	else
	{	
		KdInterior *interior = dynamic_cast<KdInterior *>(node);
	
		stream << "<Interior plane=\"" << interior->mPosition << " " 
			   << interior->mAxis << "\">" << endl;

		ExportNode(interior->mBack, stream);
		ExportNode(interior->mFront, stream);

		stream << "</Interior>" << endl;
	}
}


struct KdObjectsTraversalData
{
	KdNode *node;
	ObjectContainer *objects;
};


void OspTree::InsertObjects(KdNode *node, const ObjectContainer &objects)
{
	stack<KdObjectsTraversalData> tStack;

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

		KdNode *node = tData.node;
		
		if (node->IsLeaf())
		{
			KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);

			ObjectContainer::const_iterator oit, oit_end = tData.objects->end();

			for (oit = tData.objects->begin(); oit != oit_end; ++ oit)
			{
				leaf->mObjects.push_back(*oit);
			}
		}
		else // interior
		{
			KdObjectsTraversalData frontData, backData;
			KdInterior *interior = dynamic_cast<KdInterior *>(node);

			frontData.objects = new ObjectContainer();
			backData.objects = new ObjectContainer();

			ObjectContainer::const_iterator oit, oit_end = tData.objects->end();
			
		    for (oit = tData.objects->begin(); oit != oit_end; ++ oit) 
			{
				Intersectable *object = *oit;
		
				// determine the side of this ray with respect to the plane
				const AxisAlignedBox3 box = object->GetBox();

				if (box.Max(interior->mAxis) >= interior->mPosition)
				{
					frontData.objects->push_back(object);
				}

				if (box.Min(interior->mAxis) < interior->mPosition)
				{
					backData.objects->push_back(object);
				}
			}

			tStack.push(backData);
			tStack.push(frontData);
		}

		DEL_PTR(tData.objects);
	}
}


SubdivisionCandidate * OspTree::PrepareConstruction(const VssRayContainer &sampleRays,
													const ObjectContainer &objects,
													AxisAlignedBox3 *forcedObjectSpace,
													RayInfoContainer &rays)
{
	// store pointer to this tree
	OspSubdivisionCandidate::sOspTree = this;
	mOspStats.nodes = 1;
	
	// compute bounding box from objects
	ComputeBoundingBox(objects, forcedObjectSpace);

	mTermMinProbability *= mBoundingBox.GetVolume();
	mGlobalCostMisses = 0;

	//-- create new root
	KdLeaf *kdleaf = new KdLeaf(NULL, 0);
	kdleaf->mObjects = objects;
	mRoot = kdleaf;
	
	// get clipped rays
	PreprocessRays(kdleaf, sampleRays, rays);

	// probabilty is voume of all "seen" view cells
#if 1
	const float prop = EvalViewCellsVolume(kdleaf, rays);
#else
	const float prop = GetBoundingBox().GetVolume();
#endif

	//-- add first candidate for object space partition

	// create osp traversal data
	OspTraversalData oData(kdleaf, 0, &rays, 0, prop, mBoundingBox);

	//-- first split candidate
	OspSubdivisionCandidate *oSubdivisionCandidate = 
		new OspSubdivisionCandidate(oData);

	UpdateViewCellsPvs(kdleaf, rays);

	EvalSubdivisionCandidate(*oSubdivisionCandidate);

	mTotalCost = (float)objects.size() * prop / 
		mVspTree->GetBoundingBox().GetVolume();

	EvalSubdivisionStats(*oSubdivisionCandidate);

	return oSubdivisionCandidate;
}


bool OspTree::AddLeafToPvs(KdLeaf *leaf, 
						   ViewCell *vc, 
						   const float pdf, 
						   float &contribution)
{
	// add kd intersecable to pvs
	KdIntersectable *kdObj = GetOrCreateKdIntersectable(leaf);
	
	return vc->AddPvsSample(kdObj, pdf, contribution);
}


}