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


#include "ViewCell.h"
#include "Mesh.h"
#include "Intersectable.h"
#include "KdTree.h"
#include "Triangle3.h"
#include "common.h"
#include "Environment.h"
#include "ViewCellsManager.h"
#include "Exporter.h"
#include "BvHierarchy.h"



namespace GtpVisibilityPreprocessor {



template <typename T> class myless
{
public:
	bool operator() (T v1, T v2) const
	{
		return (v1->GetMergeCost() < v2->GetMergeCost());
	}
};


typedef priority_queue<ViewCell *, vector<ViewCell *>, 
					   myless<vector<ViewCell *>::value_type> > TraversalQueue;

  int ViewCell::sMailId = 10000;//2147483647;
  int ViewCell::sReservedMailboxes = 1;


float MergeCandidate::sRenderCostWeight = 0;


// pvs penalty can be different from pvs size
inline static float EvalPvsPenalty(const float pvs, 
								   const float lower, 
								   const float upper)
{
	// clamp to minmax values
	if (pvs < lower)
		return (float)lower;
	if (pvs > upper)
		return (float)upper;

	return (float)pvs;
}

/** Counts contribution of the view cell to the pvs.
*/
inline int CountPvsContribution(ViewCell *vc)
{
	int count = 0;

	ObjectPvsIterator pit = vc->GetPvs().GetIterator();

	while (pit.HasMoreEntries())
	{
		Intersectable *obj = pit.Next();

		if (!obj->Mailed())
		{
			obj->Mail();
			++ count;
		}
	}

	return count;
}


/// Fast computation of merged pvs size
static float ComputeMergedPvsCost(const ObjectPvs &pvs1, 
								  const ObjectPvs &pvs2)
{
	// add first pvs
	float pvs = (float)pvs1.GetSize();

	Intersectable::NewMail();

	// mail all objects in first pvs
	ObjectPvsIterator pit = pvs1.GetIterator();

	while (pit.HasMoreEntries())
	{
		pit.Next()->Mail();
	}

	pit = pvs2.GetIterator();

	while (pit.HasMoreEntries())
	{
		if (!pit.Next()->Mailed())
			++ pvs;
	}

	return pvs;
}


ViewCell::ViewCell(): 
MeshInstance(NULL), 
mArea(-1),
mVolume(-1),
mValid(true),
mParent(NULL),
mMergeCost(0),
mPvsCost(0),
mEntriesInPvs(0),
mPvsSizeValid(false),
mFilteredPvsSize(0)
{
	mId = -100;
}

ViewCell::ViewCell(Mesh *mesh): 
MeshInstance(mesh), 
mArea(-1),
mVolume(-1),
mValid(true),
mParent(NULL),
mMergeCost(0),
mPvsCost(0),
mPvsSizeValid(false),
mFilteredPvsSize(0)
//mMailbox(0)
{
	mId = -100;
}


ViewCell::~ViewCell()
{
}


const ObjectPvs &ViewCell::GetPvs() const
{
  return mPvs;
}


ObjectPvs &ViewCell::GetPvs()
{
  return mPvs;
}


void ViewCell::SetPvs(const ObjectPvs &pvs)
{
	mPvs = pvs;
}


int ViewCell::Type() const 
{ 
	return VIEW_CELL; 
}


float ViewCell::GetVolume() const
{
	return mVolume;
}


void ViewCell::SetVolume(float volume) 
{
	mVolume = volume;
}


void ViewCell::SetMesh(Mesh *mesh)
{
	mMesh = mesh;
}


float ViewCell::GetArea() const
{
	return mArea;
}


void ViewCell::SetArea(float area) 
{
	mArea = area;
}


void ViewCell::SetColor(const RgbColor &color) 
{
	mColor = color;
}


RgbColor ViewCell::GetColor() const 
{
	return mColor;
}


void ViewCell::SetValid(const bool valid)
{
	mValid = valid;
}


bool ViewCell::GetValid() const
{
	return mValid;
}


void ViewCell::SetParent(ViewCellInterior *parent)
{
	mParent = parent;
}


bool ViewCell::IsRoot() const
{
	return !mParent;
}


ViewCellInterior *ViewCell::GetParent() const
{
	return mParent;
}


void ViewCell::SetMergeCost(const float mergeCost)
{
	mMergeCost = mergeCost;
}


float ViewCell::GetRenderCost() const
{
	return mPvsCost * GetVolume();
}


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


bool ViewCell::AddPvsSample(Intersectable *sample,
							const float pdf, 
							float &contribution)
{ 
	const bool result = mPvs.AddSample(sample, pdf);
	// have to recompute pvs size
	mPvsSizeValid = false;

	return result;
}



/************************************************************************/
/*                class ViewCellInterior implementation                 */
/************************************************************************/


ViewCellInterior::ViewCellInterior()
{
}


ViewCellInterior::~ViewCellInterior()
{
	ViewCellContainer::const_iterator it, it_end = mChildren.end();

	for (it = mChildren.begin(); it != it_end; ++ it)
	{
		delete (*it);
	}
}


ViewCellInterior::ViewCellInterior(Mesh *mesh): 
ViewCell(mesh)
{
}


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


void ViewCellInterior::SetupChildLink(ViewCell *vc)
{
    mChildren.push_back(vc);
    vc->SetParent(this);
}


void ViewCellInterior::RemoveChildLink(ViewCell *l)
{
	// erase leaf from old view cell
	ViewCellContainer::iterator it = mChildren.begin();

	for (; (*it) != l; ++ it);
	if (it == mChildren.end())
		Debug << "error" << endl;
	else
		mChildren.erase(it);
}


void ViewCellInterior::ReplaceChildLink(ViewCell *prev, ViewCell *cur)
{
	// erase leaf from old view cell
	ViewCellContainer::iterator it = mChildren.begin();

	for (; (*it) != prev; ++ it);
	if (it == mChildren.end())
	{
		Debug << "error: child link not found" << endl;
	}
	else
	{
		(*it) = cur;
	}
}



/************************************************************************/
/*                class ViewCellsStatistics implementation              */
/************************************************************************/


void ViewCellsStatistics::Print(ostream &app) const
{
	app << "=========== View Cells Statistics ===============\n";

	app << setprecision(4);

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

	app << "#N_OVERALLPVS ( cost of the PVS )\n" << pvsCost << endl;

	//app << "#N_PVSENTRIES ( entries in the PVS)\n" << pvsEntries << endl;

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

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

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

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

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

	app << "#N_AVGLEAVES (average number of leaves per view cell )\n" << AvgLeaves() << endl;

	app << "#N_MAXLEAVES ( maximal number of leaves per view cell )\n" << maxLeaves << endl;
	
	app << "#N_INVALID ( number of invalid view cells )\n" << invalid << endl;

	app << "========== End of View Cells Statistics ==========\n";
}


/*************************************************************************/
/*                    class ViewCellsTree implementation                 */
/*************************************************************************/


ViewCellsTree::ViewCellsTree(ViewCellsManager *vcm):
mRoot(NULL),
mUseAreaForPvs(false),
mViewCellsManager(vcm),
#if 0
mViewCellsStorage(PVS_IN_INTERIORS)
#else
mViewCellsStorage(PVS_IN_LEAVES)
#endif
{
	ReadEnvironment();
	MergeCandidate::sRenderCostWeight = mRenderCostWeight;
}


ViewCellsTree::ViewCellsTree():
mRoot(NULL),
mUseAreaForPvs(false),
mViewCellsManager(NULL),
#if 0
mViewCellsStorage(PVS_IN_INTERIORS)
#else
mViewCellsStorage(PVS_IN_LEAVES)
#endif
{
	ReadEnvironment();
	MergeCandidate::sRenderCostWeight = mRenderCostWeight;
}


void ViewCellsTree::ReadEnvironment()
{
	Environment::GetSingleton()->GetBoolValue("ViewCells.Visualization.exportMergedViewCells", mExportMergedViewCells);
	Environment::GetSingleton()->GetFloatValue("ViewCells.maxStaticMemory", mMaxMemory);

	//-- merge options
	Environment::GetSingleton()->GetFloatValue("ViewCells.PostProcess.renderCostWeight", mRenderCostWeight);
	Environment::GetSingleton()->GetIntValue("ViewCells.PostProcess.minViewCells", mMergeMinViewCells);
	Environment::GetSingleton()->GetFloatValue("ViewCells.PostProcess.maxCostRatio", mMergeMaxCostRatio);
	Environment::GetSingleton()->GetBoolValue("ViewCells.PostProcess.refine", mRefineViewCells);	

	Environment::GetSingleton()->GetIntValue("ViewCells.PostProcess.maxMergesPerPass", mMaxMergesPerPass);
	Environment::GetSingleton()->GetFloatValue("ViewCells.PostProcess.avgCostMaxDeviation", mAvgCostMaxDeviation);

	Debug << "============= view cell tree options ================\n";
	Debug << "minimum view cells: " << mMergeMinViewCells << endl;
	Debug << "max cost ratio: " << mMergeMaxCostRatio << endl;
	Debug << "max memory: " << mMaxMemory << endl;
	Debug << "refining view cells: " << mRefineViewCells << endl;
	Debug << "=========== end view cell tree options ===============\n";
}


// return memory usage in MB
float ViewCellsTree::GetMemUsage() const
{
	// TODO
	return 0;
		/*(sizeof(ViewCellsTree) +
		 mBspStats.Leaves() * sizeof(BspLeaf) +
		 mBspStats.Interior() * sizeof(BspInterior) +
		 mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);*/
}


int ViewCellsTree::GetNumInitialViewCells(ViewCell *vc) const
{
	int vcSize = 0;

	stack<ViewCell *> tstack;

	tstack.push(vc);

	while (!tstack.empty())
	{
		ViewCell *vc = tstack.top();
		tstack.pop();

		if (vc->IsLeaf())
		{
			++ vcSize;
		}
		else
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
			
			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}
			
		}
	}

	return vcSize;
}


void ViewCellsTree::CollectLeaves(ViewCell *vc, ViewCellContainer &leaves) const
{
	stack<ViewCell *> tstack;

	tstack.push(vc);

	while (!tstack.empty())
	{
		ViewCell *vc = tstack.top();
		tstack.pop();

		if (vc->IsLeaf())
		{
			leaves.push_back(vc);
		}
		else
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}
		}
	}
}


ViewCellsTree::~ViewCellsTree()
{
	DEL_PTR(mRoot);
}


int ViewCellsTree::ConstructMergeTree(const VssRayContainer &rays, 
									  const ObjectContainer &objects)
{
	mNumActiveViewCells = (int)mViewCellsManager->GetViewCells().size();

	float variance = 0;
	float totalPvs = 0;
	float totalRenderCost = 0;

	//////////////////
	//-- compute statistics values of initial view cells

	mViewCellsManager->EvaluateRenderStatistics(totalRenderCost,
												mExpectedCost,
												mDeviation,
												variance,
												totalPvs,
												mAvgRenderCost);

	/////////////
	//-- fill merge queue
	vector<MergeCandidate> candidates;

	mViewCellsManager->CollectMergeCandidates(rays, candidates);

	while(!candidates.empty())
	{
		MergeCandidate mc = candidates.back();
		candidates.pop_back();
		EvalMergeCost(mc);
		mMergeQueue.push(mc);
	}

	Debug << "************************* merge ***********************************" << endl;  
	Debug << "deviation: " << mDeviation << endl;
	Debug << "avg render cost: " << mAvgRenderCost << endl;
	Debug << "expected cost: " << mExpectedCost << endl;


	ViewCellsManager::PvsStatistics pvsStats;
	mViewCellsManager->GetPvsStatistics(pvsStats);

	//static float expectedValue = pvsStats.avgPvs;
	
	//-- the current view cells are kept in this container
	//-- we start with the current view cells from the view cell manager. 
	//-- The active view cells will change with subsequent merges
	
	// todo: should rather take initial view cells
    ViewCellContainer &activeViewCells = mViewCellsManager->GetViewCells();
	
	
	ViewCell::NewMail();

	MergeStatistics mergeStats;
	mergeStats.Start();
	
	long startTime = GetTime();

	mergeStats.collectTime = TimeDiff(startTime, GetTime());
	mergeStats.candidates = (int)mMergeQueue.size();
	startTime = GetTime();

	// frequency stats are updated
	const int statsOut = 500;
	
	// passes are needed for statistics, because we don't want to record
	// every merge
	int pass = 0;
	int mergedPerPass = 0;
	float realExpectedCost = mExpectedCost;
	float realAvgRenderCost = mAvgRenderCost;
	int realNumActiveViewCells = mNumActiveViewCells;
	
	// maximal ratio of old expected render cost to expected render
	// when the the render queue has to be reset.
	int numMergedViewCells = 0;
		

	cout << "actual merge starts now ... " << endl;

	//-- use priority queue to merge leaf pairs

	while (!mMergeQueue.empty())
	{
		//-- reset merge queue if the ratio of current expected cost / real expected cost
		//   too small or after a given number of merges
		if ((mergedPerPass > mMaxMergesPerPass) ||
			(mAvgCostMaxDeviation > mAvgRenderCost / realAvgRenderCost))
		{
			Debug << "************ reset queue *****************\n"
				  << "ratios: " << mAvgCostMaxDeviation 
				  << " real avg render cost " << realAvgRenderCost << " average render cost " << mAvgRenderCost
				  << " merged per pass : " << mergedPerPass << " of maximal " << mMaxMergesPerPass << endl;

			Debug << "Values before reset: "  
				  << " erc: " << mExpectedCost 
				  << " avgrc: " << mAvgRenderCost
				  << " dev: " << mDeviation << endl;
	
			// adjust render cost
			++ pass;

			mergedPerPass = 0;
			mExpectedCost = realExpectedCost;
			mAvgRenderCost = realAvgRenderCost;
			mNumActiveViewCells = realNumActiveViewCells;
			
			const int numMergedViewCells = UpdateActiveViewCells(activeViewCells);
		
			/////////////////
			//-- reset / refine the view cells
			//-- priorities are recomputed
			//-- the candidates are put back into merge queue

			if (mRefineViewCells)
				RefineViewCells(rays, objects);
			else
				ResetMergeQueue();

			Debug << "Values after reset: "  
				  << " erc: " << mExpectedCost 
				  << " avg: " << mAvgRenderCost 
				  << " dev: " << mDeviation << endl;

			if (mExportMergedViewCells)
			{
				ExportMergedViewCells(activeViewCells, objects, numMergedViewCells);
			}
		}


		MergeCandidate mc = mMergeQueue.top();
		mMergeQueue.pop();
	
		// both view cells equal because of previous merges
		// NOTE: do I really still need this? probably cannot happen!!
		if (mc.mLeftViewCell == mc.mRightViewCell)
			continue;

		if (mc.IsValid())
		{
			ViewCell::NewMail();

			//-- update statistical values
			-- realNumActiveViewCells;
			++ mergeStats.merged;
			++ mergedPerPass;

			const float renderCostIncr = mc.GetRenderCost();
			const float mergeCostIncr = mc.GetMergeCost();

			totalRenderCost += renderCostIncr;
			mDeviation += mc.GetDeviationIncr();

						
			//-- merge the view cells of leaf1 and leaf2
			float pvsDiff;
			ViewCellInterior *mergedVc = 
				MergeViewCells(mc.mLeftViewCell, mc.mRightViewCell, pvsDiff);	

			// total render cost and deviation has changed
			// real expected cost will be larger than expected cost used for the
			// cost heuristics, but cannot recompute costs on each increase of the 
			// expected cost
			totalPvs += pvsDiff;
			realExpectedCost = totalRenderCost / (float)realNumActiveViewCells;
			realAvgRenderCost = (float)totalPvs / (float)realNumActiveViewCells;
	
			// set merge cost to this node for priority traversal
			mergedVc->SetMergeCost(totalRenderCost);
			
			// check if "siblings (back and front node of the same parent)
			if (0) ++ mergeStats.siblings;

			// set the cost for rendering a view cell
			mergedVc->SetCost(realExpectedCost);

			if ((mergeStats.merged % statsOut) == 0)
				cout << "merged " << mergeStats.merged << " view cells" << endl;

		}
		else
		{ 
			// merge candidate not valid, because one of the leaves was already
			// merged with another one => validate and reinsert into queue
			if (ValidateMergeCandidate(mc))
			{
				EvalMergeCost(mc);
				mMergeQueue.push(mc);
			}
		}
		
	}

	// adjust stats and reset queue one final time
	mExpectedCost = realExpectedCost;
	mAvgRenderCost = realAvgRenderCost;
	mNumActiveViewCells = realNumActiveViewCells;

	UpdateActiveViewCells(activeViewCells);

	// refine view cells and reset costs
	if (mRefineViewCells)
		RefineViewCells(rays, objects);
	else
		ResetMergeQueue();


	// create a root node if the merge was not done till root level,
	// else take the single node as new root
	if ((int)activeViewCells.size() > 1)
	{
		Debug << "creating root of view cell hierarchy for " 
			  << (int)activeViewCells.size() << " view cells" << endl;
		
		ViewCellInterior *root = mViewCellsManager->MergeViewCells(activeViewCells);
	
		ViewCellContainer::const_iterator it, it_end = root->mChildren.end();

		for (it = root->mChildren.begin(); it != it_end; ++ it)
			(*it)->SetParent(root);
	
		root->SetMergeCost(totalRenderCost);
		// $$JB keep this 0 temporarilly
		root->SetCost(0.0f);

		mRoot = root;
	}
	// normal case
	else if (!activeViewCells.empty())
	{
		Debug << "setting root of the merge history" << endl;
		mRoot = activeViewCells[0];
	}
	else
	{
		Debug << "big error, root is NULL" << endl;
	}
	
	//-- empty merge queue just in case
	while (!mMergeQueue.empty())
	{
		mMergeQueue.pop();
	}

	// test if computed volumes are correct
	Debug << "volume of the root view cell: " << mRoot->GetVolume() 
		  << " " << mViewCellsManager->GetViewSpaceBox().GetVolume() << endl;
	
	// TODO: delete because makes no sense here
	mergeStats.expectedRenderCost = realExpectedCost;
	mergeStats.deviation = mDeviation;

	// we want to optimize this heuristics
	mergeStats.heuristics = 
		mDeviation * (1.0f - mRenderCostWeight) + 
		mExpectedCost * mRenderCostWeight;

	mergeStats.mergeTime = TimeDiff(startTime, GetTime());
	mergeStats.Stop();
	Debug << mergeStats << endl << endl;

	// assign colors for the view cells so that at least one is always consistent
	AssignRandomColors();

	//TODO: should return sample contributions?
	return mergeStats.merged;
}


ViewCell *ViewCellsTree::GetRoot() const
{
	return mRoot;
}


void ViewCellsTree::ResetMergeQueue()
{
	cout << "reset merge queue ... ";
	
	vector<MergeCandidate> buf;
	buf.reserve(mMergeQueue.size());
			
	
	// store merge candidates in intermediate buffer
	while (!mMergeQueue.empty())
	{
		MergeCandidate mc = mMergeQueue.top();
		mMergeQueue.pop();
		
		// recalculate cost
		if (ValidateMergeCandidate(mc))
		{
			EvalMergeCost(mc);
			buf.push_back(mc);				
		}
	}

	vector<MergeCandidate>::const_iterator bit, bit_end = buf.end();

	// reinsert back into queue
	for (bit = buf.begin(); bit != bit_end; ++ bit)
	{      
		mMergeQueue.push(*bit);	
	}

	cout << "finished" << endl;
}


float ViewCellsTree::ComputeMergedPvsCost(const ObjectPvs &pvs1, 
										  const ObjectPvs &pvs2) const
{
	// computes render cost of merge
	float renderCost = 0;

	// compute new pvs size
	Intersectable::NewMail();

	ObjectPvsIterator pit = pvs1.GetIterator();

	// first mail all objects in first pvs
	while (pit.HasMoreEntries())
	{
		Intersectable *obj = pit.Next();

		obj->Mail();
		renderCost += mViewCellsManager->EvalRenderCost(obj);
	}

	ObjectPvsIterator pit2 = pvs2.GetIterator();

	while (pit2.HasMoreEntries())
	{
		Intersectable *obj = pit2.Next();
		
		// test if object already considered	
		if (!obj->Mailed())
		{
			renderCost += mViewCellsManager->EvalRenderCost(obj);
		}
	}

	return renderCost;
}


int ViewCellsTree::UpdateActiveViewCells(ViewCellContainer &viewCells)
{
	int numMergedViewCells = 0;

	Debug << "updating active vc: " << (int)viewCells.size() << endl;
	
	// find all already merged view cells and remove them from the
	// container view cells
		
	// sort out all view cells which are not active anymore, i.e., they
	// were already part of a merge
	int i = 0;

	ViewCell::NewMail();

	while (1)
	{
		// remove all merged view cells from end of the vector
		while (!viewCells.empty() && (viewCells.back()->GetParent()))
		{
			viewCells.pop_back();
		}

		// all merged view cells have been found
		if (i >= (int)viewCells.size()) 
			break;

		// already merged this view cell, put it to end of vector
		if (viewCells[i]->GetParent())
			swap(viewCells[i], viewCells.back());
		
		// mail view cell that it has not been merged
		viewCells[i]->Mail();

		// increase loop counter
		++ i;
	}


	// add new view cells to container only if they don't have been 
	// merged in the mean time
	ViewCellContainer::const_iterator ait, ait_end = mMergedViewCells.end();

	for (ait = mMergedViewCells.begin(); ait != ait_end; ++ ait)
	{
		ViewCell *vc = mMergedViewCells.back();
		if (!vc->GetParent() && !vc->Mailed())
		{
			vc->Mail();
			viewCells.push_back(vc);
			++ numMergedViewCells;
		}
	}

	// dispose old merged view cells
	mMergedViewCells.clear();

	// update standard deviation
	ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
	
	mDeviation = 0;

	for (vit = viewCells.begin(); vit != vit_end; ++ vit)
	{
		const float lower = (float)mViewCellsManager->GetMinPvsSize();
		const float upper = (float)mViewCellsManager->GetMaxPvsSize();

		const float penalty = EvalPvsPenalty((*vit)->GetPvs().EvalPvsCost(), lower, upper);
		
		mDeviation += fabs(mAvgRenderCost - penalty);
	}

	mDeviation /= (float)viewCells.size();
	
	return numMergedViewCells;
}


void ViewCellsTree::ExportMergedViewCells(ViewCellContainer &viewCells, 
										  const ObjectContainer &objects,
										  const int numMergedViewCells)
{
	

	char s[64];

	sprintf(s, "merged_viewcells%07d.x3d", (int)viewCells.size());
	Exporter *exporter = Exporter::GetExporter(s);

	if (exporter)
	{
		cout << "exporting " << (int)viewCells.size() << " merged view cells ... ";
		exporter->ExportGeometry(objects);
		//Debug << "vc size " << (int)viewCells.size() << " merge queue size: " << (int)mMergeQueue.size() << endl;
		ViewCellContainer::const_iterator it, it_end = viewCells.end();

		int i = 0;
		for (it = viewCells.begin(); it != it_end; ++ it)
		{
			Material m;
			// assign special material to new view cells
			// new view cells are on the back of container
			if (i ++ >= (viewCells.size() - numMergedViewCells))
			{
				//m = RandomMaterial();
				m.mDiffuseColor.r = RandomValue(0.5f, 1.0f);
				m.mDiffuseColor.g = RandomValue(0.5f, 1.0f);
				m.mDiffuseColor.b = RandomValue(0.5f, 1.0f);
			}
			else
			{
				float col = RandomValue(0.1f, 0.4f);
				m.mDiffuseColor.r = col;
				m.mDiffuseColor.g = col;
				m.mDiffuseColor.b = col;
			}

			exporter->SetForcedMaterial(m);
			mViewCellsManager->ExportViewCellGeometry(exporter, *it, NULL, NULL);
		}

		delete exporter;
		cout << "finished" << endl;
	}
}


// TODO: should be done in view cells manager
ViewCellInterior *ViewCellsTree::MergeViewCells(ViewCell *left, 
												ViewCell *right, 
												float &pvsDiff) //const
{
	// create merged view cell
	ViewCellInterior *vc = 
		mViewCellsManager->MergeViewCells(left, right);

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

	// set to the new parent view cell
	left->SetParent(vc);
	right->SetParent(vc);

	
	if (mUseAreaForPvs)
	{
		// set new area of view cell
		// not not correct, but costly to compute real area!!
		vc->SetArea(left->GetArea() + right->GetArea());
	}
	else
	{	// set new volume of view cell
		vc->SetVolume(left->GetVolume() + right->GetVolume());
	}

	
	// important so other merge candidates sharing this view cell
	// are notified that the merge cost must be updated!!
	vc->Mail();

	const float pvs1 = left->GetPvs().EvalPvsCost();
	const float pvs2 = right->GetPvs().EvalPvsCost();


	// the new view cells are stored in this container
	mMergedViewCells.push_back(vc);

	pvsDiff = vc->GetPvs().EvalPvsCost() - pvs1 - pvs2;


	// don't store pvs in interior cells, just a scalar
	if (mViewCellsStorage == PVS_IN_LEAVES)
	{
		// set scalars
		mViewCellsManager->UpdateScalarPvsSize(left, 
				left->GetPvs().EvalPvsCost(), 
				left->GetPvs().GetSize());
			
		// remove pvs, we don't store interior pvss
		if (!left->IsLeaf())
		{
			left->GetPvs().Clear();
		}

		// set scalars
		mViewCellsManager->UpdateScalarPvsSize(right, 
			right->GetPvs().EvalPvsCost(), 
			right->GetPvs().GetSize());

		right->mPvsSizeValid = true;
		
		// remove pvs, we don't store interior pvss
		if (!right->IsLeaf())
		{
			right->GetPvs().Clear();
		}
	}

	return vc;
}


int ViewCellsTree::RefineViewCells(const VssRayContainer &rays, 
								   const ObjectContainer &objects)
{
	Debug << "refining " << (int)mMergeQueue.size() << " candidates " << endl;

	// intermediate buffer for shuffled view cells
	vector<MergeCandidate> buf;
	buf.reserve(mMergeQueue.size());
			
	// Use priority queue of remaining leaf pairs 
	// The candidates either share the same view cells or
	// are border leaves which share a boundary.
	// We test if they can be shuffled, i.e.,
	// either one leaf is made part of one view cell or the other
	// leaf is made part of the other view cell. It is tested if the
	// remaining view cells are "better" than the old ones.
	
	const int numPasses = 3;
	int pass = 0;
	int passShuffled = 0;
	int shuffled = 0;
	int shuffledViewCells = 0;

	ViewCell::NewMail();
	
	while (!mMergeQueue.empty())
	{
		MergeCandidate mc = mMergeQueue.top();
		mMergeQueue.pop();

		// both view cells equal or already shuffled
		if ((mc.GetLeftViewCell() == mc.GetRightViewCell()) ||
			mc.GetLeftViewCell()->IsLeaf() || mc.GetRightViewCell()->IsLeaf())
		{			
			continue;
		}

		// candidate for shuffling
		const bool wasShuffled = ShuffleLeaves(mc);
		
		// shuffled or put into other queue for further refine
		if (wasShuffled)
		{
			++ passShuffled;

			if (!mc.GetLeftViewCell()->Mailed())
			{
				mc.GetLeftViewCell()->Mail();
				++ shuffledViewCells;
			}
			if (!mc.GetRightViewCell()->Mailed())
			{
				mc.GetRightViewCell()->Mail();
				++ shuffledViewCells;
			}
		}

		// put back into intermediate vector
		buf.push_back(mc);
	}


	//-- in the end, the candidates must be in the mergequeue again
	//   with the correct cost

	cout << "reset merge queue ... ";
	
	
	vector<MergeCandidate>::const_iterator bit, bit_end = buf.end();
	
	for (bit = buf.begin(); bit != bit_end; ++ bit)
	{    
		MergeCandidate mc = *bit;
		// recalculate cost
		if (ValidateMergeCandidate(mc))
		{
			EvalMergeCost(mc);
			mMergeQueue.push(mc);	
		}
	}

	cout << "finished" << endl;

	return shuffledViewCells;
}


inline int AddedPvsSize(const ObjectPvs &pvs1, const ObjectPvs &pvs2)
{
	ObjectPvs interPvs;
	ObjectPvs::Merge(interPvs, pvs1, pvs2);

	return (int)interPvs.GetSize();
}


// recomputes pvs size minus pvs of leaf l
#if 0
inline int SubtractedPvsSize(BspViewCell *vc, BspLeaf *l, const ObjectPvs &pvs2)
{
	ObjectPvs pvs;
	vector<BspLeaf *>::const_iterator it, it_end = vc->mLeaves.end();
	for (it = vc->mLeaves.begin(); it != vc->mLeaves.end(); ++ it)
		if (*it != l)
			pvs.AddPvs(*(*it)->mPvs);
	return pvs.GetSize();
}
#endif


// computes pvs1 minus pvs2
inline int SubtractedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
{
	return pvs1.SubtractPvs(pvs2);
}


float ViewCellsTree::EvalShuffleCost(ViewCell *leaf, 
									 ViewCellInterior *vc1, 
									 ViewCellInterior *vc2) const
{
	//const int pvs1 = SubtractedPvsSize(vc1, leaf, *leaf->mPvs);
	const float pvs1 = (float)SubtractedPvsSize(vc1->GetPvs(), leaf->GetPvs());
	const float pvs2 = (float)AddedPvsSize(vc2->GetPvs(), leaf->GetPvs());

	const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
	const float upperPvsLimit = (float)mViewCellsManager->GetMaxPvsSize();

	const float pvsPenalty1 = 
		EvalPvsPenalty(pvs1, lowerPvsLimit, upperPvsLimit);

	const float pvsPenalty2 =
		EvalPvsPenalty(pvs2, lowerPvsLimit, upperPvsLimit);


	// don't shuffle leaves with pvs > max
	if (0 && (pvs1 + pvs2 > mViewCellsManager->GetMaxPvsSize()))
	{
		return 1e20f;
	}

	float p1, p2;

    if (mUseAreaForPvs)
	{
		p1 = vc1->GetArea() - leaf->GetArea();
		p2 = vc2->GetArea() + leaf->GetArea();
	}
	else
	{
		p1 = vc1->GetVolume() - leaf->GetVolume();
		p2 = vc2->GetVolume() + leaf->GetVolume();
	}

	const float renderCost1 = pvsPenalty1 * p1;
	const float renderCost2 = pvsPenalty2 * p2;

	float dev1, dev2;

	if (1)
	{
		dev1 = fabs(mAvgRenderCost - pvsPenalty1);
		dev2 = fabs(mAvgRenderCost - pvsPenalty2);
	}
	else
	{
		dev1 = fabs(mExpectedCost - renderCost1);
		dev2 = fabs(mExpectedCost - renderCost2);
	}
	
	return mRenderCostWeight * (renderCost1 + renderCost2) +
		  (1.0f - mRenderCostWeight) * (dev1 + dev2) / (float)mNumActiveViewCells;
}


void ViewCellsTree::ShuffleLeaf(ViewCell *leaf,
								ViewCellInterior *vc1, 
								ViewCellInterior *vc2) const
{
	// compute new pvs and area
	// TODO change
	vc1->GetPvs().SubtractPvs(leaf->GetPvs());
	
	ObjectPvs interPvs;
	ObjectPvs::Merge(interPvs, vc2->GetPvs(), leaf->GetPvs());
	vc2->SetPvs(interPvs);

	if (mUseAreaForPvs)
	{
		vc1->SetArea(vc1->GetArea() - leaf->GetArea());
		vc2->SetArea(vc2->GetArea() + leaf->GetArea());
	}
	else
	{
		vc1->SetVolume(vc1->GetVolume() - leaf->GetVolume());
		vc2->SetVolume(vc2->GetVolume() + leaf->GetVolume());
	}

	
	ViewCellInterior *p = static_cast<ViewCellInterior *>(leaf->GetParent());

	p->RemoveChildLink(leaf);
	vc2->SetupChildLink(leaf);
}


bool ViewCellsTree::ShuffleLeaves(MergeCandidate &mc) const
{
	float cost1, cost2;

	ViewCellInterior *vc1 = static_cast<ViewCellInterior *>(mc.GetLeftViewCell());
	ViewCellInterior *vc2 = static_cast<ViewCellInterior *>(mc.GetRightViewCell());

	ViewCell *leaf1 = mc.GetInitialLeftViewCell();
	ViewCell *leaf2 = mc.GetInitialRightViewCell();

	//-- first test if shuffling would decrease cost
	cost1 = GetCostHeuristics(vc1);
	cost2 = GetCostHeuristics(vc2);

	const float oldCost = cost1 + cost2;
	
	float shuffledCost1 = Limits::Infinity;
	float shuffledCost2 = Limits::Infinity;

	if (leaf1)
		shuffledCost1 = EvalShuffleCost(leaf1, vc1, vc2);	
	if (leaf2)
		shuffledCost2 = EvalShuffleCost(leaf2, vc2, vc1);

	// if cost of shuffle is less than old cost => shuffle
	if ((oldCost <= shuffledCost1) && (oldCost <= shuffledCost2))
		return false;
	
	if (shuffledCost1 < shuffledCost2)
	{
		if (leaf1)
			ShuffleLeaf(leaf1, vc1, vc2);
		mc.mInitialLeftViewCell = NULL;
	}
	else
	{
		if (leaf2)
			ShuffleLeaf(leaf2, vc2, vc1);
		mc.mInitialRightViewCell = NULL;
	}

	return true;
}


float ViewCellsTree::GetVariance(ViewCell *vc) const
{
	const float upper = (float)mViewCellsManager->GetMaxPvsSize();
	const float lower = (float)mViewCellsManager->GetMinPvsSize();

	if (1)
	{
		const float penalty = 
			EvalPvsPenalty(GetPvsCost(vc), lower, upper);

		return (mAvgRenderCost - penalty) * (mAvgRenderCost - penalty) / 
			(float)mNumActiveViewCells;
	}

    const float leafCost = GetRenderCost(vc);
	return (mExpectedCost - leafCost) * (mExpectedCost - leafCost);
}


float ViewCellsTree::GetDeviation(ViewCell *vc) const
{
	const float upper = (float)mViewCellsManager->GetMaxPvsSize();
	const float lower = (float)mViewCellsManager->GetMinPvsSize();

	if (1)
	{
		const float penalty = EvalPvsPenalty(GetPvsCost(vc), lower, upper);
		return fabs(mAvgRenderCost - penalty) / (float)mNumActiveViewCells;
	}

    const float renderCost = GetRenderCost(vc);
	return fabs(mExpectedCost - renderCost);
}


float ViewCellsTree::GetRenderCost(ViewCell *vc) const
{
	if (mUseAreaForPvs)
	{
		return vc->GetPvs().EvalPvsCost() * vc->GetArea();
	}

	return vc->GetPvs().EvalPvsCost() * vc->GetVolume();
}


float ViewCellsTree::GetCostHeuristics(ViewCell *vc) const
{
	return GetRenderCost(vc) * mRenderCostWeight + 
		   GetDeviation(vc) * (1.0f - mRenderCostWeight);
}


bool ViewCellsTree::ValidateMergeCandidate(MergeCandidate &mc) const
{
	// propagate up so we have only the active view cells
	while (mc.mLeftViewCell->mParent)
	{
		mc.mLeftViewCell = mc.mLeftViewCell->mParent;
	}

	while (mc.mRightViewCell->mParent)
	{
		mc.mRightViewCell = mc.mRightViewCell->mParent;
	}

	// this view cell was already merged
	//return mc.mLeftViewCell && (mc.mLeftViewCell != mc.mRightViewCell);
	return mc.mLeftViewCell != mc.mRightViewCell;
}


void ViewCellsTree::EvalMergeCost(MergeCandidate &mc) const
{
	///////////////////
	//-- compute pvs difference

	float newPvs;
	
	// not valid if not using const cost per object!!
	newPvs = ComputeMergedPvsCost(mc.mLeftViewCell->GetPvs(), mc.mRightViewCell->GetPvs());


	const float newPenalty = EvalPvsPenalty(newPvs, 
											(float)mViewCellsManager->GetMinPvsSize(),
											(float)mViewCellsManager->GetMaxPvsSize());

	ViewCell *vc1 = mc.mLeftViewCell;
	ViewCell *vc2 = mc.mRightViewCell;

	//-- compute ratio of old cost
	//-- (i.e., added size of left and right view cell times pvs size)
	//-- to new rendering cost (i.e, size of merged view cell times pvs size)
	const float oldCost = GetRenderCost(vc1) + GetRenderCost(vc2);

    const float newCost = mUseAreaForPvs ? 
		(float)newPenalty * (vc1->GetArea() + vc2->GetArea()) :
		(float)newPenalty * (vc1->GetVolume() + vc2->GetVolume());


	// strong penalty if pvs size too large
	if (0 && (newPvs > mViewCellsManager->GetMaxPvsSize()))
	{
		mc.mRenderCost = 1e20f;
	}
	else
	{
		mc.mRenderCost = (newCost - oldCost) / 
			mViewCellsManager->GetViewSpaceBox().GetVolume();
	}	
	
	///////////////////////////
	//-- merge cost also takes deviation into account

	float newDev, oldDev;

	if (1)
		newDev = fabs(mAvgRenderCost - newPenalty) / (float)mNumActiveViewCells;
	else
		newDev = fabs(mExpectedCost - newCost) / (float)mNumActiveViewCells;
	
	oldDev = GetDeviation(vc1) + GetDeviation(vc2);

	// compute deviation increase
	mc.mDeviationIncr = newDev - oldDev;
	
	//Debug << "render cost: " << mc.mRenderCost * mRenderCostWeight << endl;
	//Debug << "standard deviation: " << mc.mDeviationIncr * mRenderCostWeight << endl;
}


void ViewCellsTree::SetViewCellsStorage(int stype)
{
	if (mViewCellsStorage == stype)
		return;

	// TODO
	switch (stype)
	{
	case COMPRESSED:
		CompressViewCellsPvs(mRoot);
		break;
	default:
		break;
	}

	mViewCellsStorage = stype;
}


void ViewCellsTree::CompressViewCellsPvs(ViewCell *root)
{
	if (!root->IsLeaf())
	{
		ViewCellInterior *interior = static_cast<ViewCellInterior *>(root);

        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
		
		// compress child sets first
		for (it = interior->mChildren.begin(); it != it_end; ++ it)
		{
			CompressViewCellsPvs(*it);
		}

		/*cout << "***************\nbefore pulling up:\n";
		for (it = interior->mChildren.begin(); it != it_end; ++ it)
		{
			cout << "vc:\n" << (*it)->GetPvs() << endl;
		}*/

		// compress root node
		PullUpVisibility(interior);

		/*cout << "after pulling up:\n";
		cout << "interior:\n" << interior->GetPvs() << endl;

		for (it = interior->mChildren.begin(); it != it_end; ++ it)
		{
			cout << "vc:\n" << (*it)->GetPvs() << endl;
		}*/
	}
}


void ViewCellsTree::UpdateStats(ofstream &stats,
								const ViewCellsTreeStats &vcStats)
{
	 stats << "#Pass\n" << vcStats.mPass << endl
		   << "#Splits\n" << vcStats.mNumViewCells << endl 
		   << "#RenderCostDecrease\n" << vcStats.mRenderCostDecrease << endl // TODO
		   << "#TotalRenderCost\n" << vcStats.mTotalRenderCost << endl
		   << "#CurrentPvs\n" << vcStats.mCurrentPvsCost << endl
		   << "#ExpectedCost\n" << vcStats.mExpectedCost << endl
		   << "#AvgRenderCost\n" << vcStats.mAvgRenderCost << endl
		   << "#Deviation\n" << vcStats.mDeviation << endl
		   << "#TotalPvs\n" << vcStats.mTotalPvsCost << endl
		   << "#TotalEntriesInPvs\n" << vcStats.mEntriesInPvs << endl
		   << "#Memory\n" << vcStats.mMemoryCost << endl
		   << "#PvsSizeDecrease\n" << vcStats.mPvsSizeDecr << endl
		   << "#Volume\n" << vcStats.mVolume << endl
		   << endl;
}


void ViewCellsTree::ExportStats(const string &mergeStats)
{
	TraversalQueue tqueue;
	tqueue.push(mRoot);

	ViewCellsTreeStats vcStats;
	vcStats.Reset();

	//cout << "exporting stats ... " << endl;
	vcStats.mNumViewCells = 1;
	
	const AxisAlignedBox3 box = mViewCellsManager->GetViewSpaceBox();
	const float vol = box.GetVolume();

	const float rootPvs = GetPvsCost(mRoot);
	const int rootEntries = GetPvsEntries(mRoot);
	
	vcStats.mTotalPvsCost = rootPvs;
	vcStats.mTotalRenderCost = rootPvs;
	vcStats.mAvgRenderCost = rootPvs;
	vcStats.mExpectedCost = rootPvs;
	
	vcStats.mEntriesInPvs = rootEntries;

	ofstream stats;
	stats.open(mergeStats.c_str());

	vcStats.mMemoryCost = (float)vcStats.mEntriesInPvs * ObjectPvs::GetEntrySize();

	/////////////
	//-- first view cell

	UpdateStats(stats, vcStats);
				
		
	//-- go through tree in the order of render cost decrease
	//-- which is the same order as the view cells were merged
	//-- or the reverse order of subdivision for subdivision-only 
	//-- view cell hierarchies.

	while (!tqueue.empty())
	{
		ViewCell *vc = tqueue.top();
		tqueue.pop();

		if (!vc->IsLeaf()) 
		{	
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

			const float parentCost = GetPvsCost(interior);
			const int parentPvsEntries = GetPvsEntries(interior);
            const float parentExpCost = (float)parentCost * interior->GetVolume();

			float childExpCost = 0;
			float childCost = 0;
			int childPvsEntries = 0;

			-- vcStats.mNumViewCells;

			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				ViewCell *vc = *it;

				const float pvsCost = GetPvsCost(vc);
				const int pvsEntries = GetPvsEntries(vc);

				childExpCost += childCost * vc->GetVolume();
				childCost += pvsCost;
				childPvsEntries += pvsEntries;

				tqueue.push(vc);
				++ vcStats.mNumViewCells;
			}

			// update stats for this view cell
			const float costDecr = (parentExpCost - childExpCost) / vol;

			vcStats.mTotalRenderCost -= costDecr;
			vcStats.mTotalPvsCost += childCost - parentCost;
			vcStats.mEntriesInPvs += childPvsEntries - parentPvsEntries;

			vcStats.mExpectedCost = vcStats.mTotalRenderCost / (float)vcStats.mNumViewCells;
			vcStats.mAvgRenderCost = vcStats.mTotalPvsCost / (float)vcStats.mNumViewCells;

			vcStats.mMemoryCost = vcStats.mEntriesInPvs * ObjectPvs::GetEntrySize();

			UpdateStats(stats, vcStats);
		}
	}

	stats.close();
}


void ViewCellsTree::SetRoot(ViewCell *root)
{
	mRoot = root;
}


void ViewCellsTree::CollectBestViewCellSet(ViewCellContainer &viewCells, 
										   const int numViewCells)
{
	TraversalQueue tqueue;
	tqueue.push(mRoot);
	
	while (!tqueue.empty())
	{
		ViewCell *vc = tqueue.top();
		tqueue.pop();

		// save the view cells if it is a leaf or if enough view cells have already been traversed
		// because of the priority queue, this will be the optimal set of v
		if (vc->IsLeaf() || ((viewCells.size() + tqueue.size() + 1) >= numViewCells))
		{
			viewCells.push_back(vc);
		}
		else 
		{	
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tqueue.push(*it);
			}
		}
	}
}


void ViewCellsTree::PullUpVisibility(ViewCellInterior *interior)
{
	const int mail = (int)interior->mChildren.size();
	Intersectable::NewMail(mail);

	ViewCellContainer::const_iterator cit, cit_end = interior->mChildren.end();

    // mail all objects in the leaf sets
	// we are interested in the objects which are present in all leaves
	// => count how often an object is part of a child set
	for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
	{
		ViewCell *vc = *cit;
		ObjectPvsIterator pit = vc->GetPvs().GetIterator();

		while (pit.HasMoreEntries())
		{		
			pit.Next()->Mail();
		}
	}

	for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
	{
		ViewCell *vc = *cit;

		ObjectPvsIterator pit = vc->GetPvs().GetIterator();

		while (pit.HasMoreEntries())
		{		
			pit.Next()->IncMail();
		}
	}

	// reset pvs
    interior->GetPvs().Clear(false);
		
	PvsData pvsData;

	// only the objects which are present in all leaf pvs 
	// should remain in the parent pvs
	// these are the objects which have been mailed in all children
	for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
	{
		ViewCell *vc = *cit;
		
		ObjectPvsIterator pit = vc->GetPvs().GetIterator();

		while (pit.HasMoreEntries())
		{		
			Intersectable *obj = pit.Next(pvsData);

			if (obj->Mailed(mail))
			{	
				interior->GetPvs().AddSample(obj, pvsData.mSumPdf);
			}
		}
	}

	// delete entries which are pulled up
	// note: could be inefficent, rather gather unmailed pvs and reassign
	for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
	{
		ViewCell *vc = *cit;
		ObjectPvsIterator pit = vc->GetPvs().GetIterator();

		ObjectPvs newPvs;
		
		while (pit.HasMoreEntries())
		{		
			Intersectable *obj = pit.Next(pvsData);

			if (!obj->Mailed(mail))
			{
				newPvs.AddSampleDirty(obj, pvsData.mSumPdf);
			}
		}

		vc->SetPvs(newPvs);
	}
}


// TODO matt: implement this function for different storing methods
void ViewCellsTree::GetPvs(ViewCell *vc, ObjectPvs &pvs) const
{
	// pvs is stored in each cell => just return pvs
	if (mViewCellsStorage == PVS_IN_INTERIORS)
	{
		pvs = vc->GetPvs();
		return;
	}

	////////////////
	//-- pvs is not stored with the interiors => reconstruct
	ViewCell *root = vc;
	
	// add pvs from this view cell to root
	while (root->GetParent())
	{
		root = root->GetParent();
		pvs.MergeInPlace(root->GetPvs());
	}

	// add pvs from leaves
	stack<ViewCell *> tstack;
	tstack.push(vc);

	while (!tstack.empty())
	{
		vc = tstack.top();
		tstack.pop();
	
		// add newly found pvs to merged pvs
		pvs.MergeInPlace(vc->GetPvs());
		
		if (!vc->IsLeaf()) // interior cells: go down to leaf level
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}		
		}
	}
}


float ViewCellsTree::GetPvsCostForLeafStorage(ViewCell *vc) const
{
	// pvs is always stored directly in leaves
	if (vc->IsLeaf())
	{
		return vc->GetPvs().EvalPvsCost();
	}
	
	// interior nodes: pvs is either stored as a 
	// scalar or has to be reconstructed from the leaves
	// the stored pvs size is the valid pvs size => just return scalar
	if (vc->mPvsSizeValid)
	{
		return vc->mPvsCost;
	}
	
	// if no valid pvs size stored as a scalar => 
	// compute current pvs size of interior from it's leaf nodes
	ViewCellContainer leaves;
	CollectLeaves(vc, leaves);

	ViewCellContainer::const_iterator it, it_end = leaves.end();

	ObjectPvs newPvs;

	// sum up uniquely found intersectables
	for (it = leaves.begin(); it != it_end; ++ it)
	{
		newPvs.MergeInPlace((*it)->GetPvs());
	}

	return newPvs.EvalPvsCost();
}


int ViewCellsTree::GetEntriesInPvsForLeafStorage(ViewCell *vc) const
{
	// pvs is always stored directly in leaves
	if (vc->IsLeaf())
	{
		return vc->GetPvs().GetSize();
	}
	
	// interior node

	// interior nodes: pvs is either stored as a scalar or
	// has to be reconstructed from the leaves

	// the stored pvs size is the valid pvs size => just return scalar
	if (vc->mPvsSizeValid)
	{
		return vc->mEntriesInPvs;
	}
	
	int pvsSize = 0;

	// if no valid pvs size stored as a scalar => 
	// compute current pvs size of interior from it's leaf nodes
	ViewCellContainer leaves;
	CollectLeaves(vc, leaves);

	ViewCellContainer::const_iterator it, it_end = leaves.end();
	Intersectable::NewMail();

	// sum different intersectables
	for (it = leaves.begin(); it != it_end; ++ it)
	{
		ObjectPvsIterator oit = (*it)->GetPvs().GetIterator();

		while (oit.HasMoreEntries())
		{
			Intersectable *intersect = oit.Next();

			if (!intersect->Mailed())
			{
				intersect->Mail();
				++ pvsSize;
			}
		}
	}

	return pvsSize;
}


float ViewCellsTree::GetPvsCostForCompressedStorage(ViewCell *vc) const
{
	float pvsCost = 0;

	/////////////
	//-- compressed pvs

	// the stored pvs size is the valid pvs size => just return scalar
	if (vc->mPvsSizeValid)
	{
		pvsCost = vc->mPvsCost;
	}

	// if no pvs size stored, compute new one
	ViewCell *root = vc;
	
	// also add pvs from this view cell to root
	while (root->GetParent())
	{
		root = root->GetParent();
	
		// matt: bug! must evaluate kd pvss also
		pvsCost += CountPvsContribution(root);
	}

	stack<ViewCell *> tstack;
	tstack.push(vc);

	while (!tstack.empty())
	{
		vc = tstack.top();
		tstack.pop();

		// matt: bug! must evaluate kd pvss also
		pvsCost += CountPvsContribution(vc);

		if (!vc->IsLeaf())
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}		
		}
	}

	return pvsCost;
}


int ViewCellsTree::GetEntriesInPvsForCompressedStorage(ViewCell *vc) const
{
	int pvsSize = 0;

	/////////////////
	//-- compressed pvs

	// the stored pvs size is the valid pvs size => just return scalar
	if (vc->mPvsSizeValid)
	{
		pvsSize = vc->mEntriesInPvs;
	}

	// if no pvs size stored, compute new one
	ViewCell *root = vc;
	
	// also add pvs from this view cell to root
	while (root->GetParent())
	{
		root = root->GetParent();
		// count the pvs entries different from the already found ones	
		pvsSize += CountPvsContribution(root);
	}

	stack<ViewCell *> tstack;
	tstack.push(vc);

	while (!tstack.empty())
	{
		vc = tstack.top();
		tstack.pop();
	
		// count the pvs entries different from the already found ones	
		pvsSize += CountPvsContribution(vc);

		if (!vc->IsLeaf())
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}		
		}
	}

	return pvsSize;
}


float ViewCellsTree::GetPvsCost(ViewCell *vc) const
{
	float pvsCost = 0;
	Intersectable::NewMail();

	////////////////////////////////////////////////
	// for interiors, pvs can be stored using different methods
	//

	switch (mViewCellsStorage)
	{
	case PVS_IN_LEAVES: 
		// pvs is stored only in leaves
		pvsCost = GetPvsCostForLeafStorage(vc);			
		break;
	case COMPRESSED:
		pvsCost = GetPvsCostForCompressedStorage(vc);
		break;
	case PVS_IN_INTERIORS:
	default:
		// pvs is stored consistently in the tree up
		//to the root just return pvs size 
		pvsCost = vc->GetPvs().EvalPvsCost();	
		break;
	}

	return pvsCost;  
}


int ViewCellsTree::GetPvsEntries(ViewCell *vc) const
{
	int pvsSize = 0;

	Intersectable::NewMail();

	////////////////////////////////////////////////
	// for interiors, pvs can be stored using different methods
	
	switch (mViewCellsStorage)
	{
	case PVS_IN_LEAVES:
		{
			//-- store pvs only in leaves
			pvsSize = GetEntriesInPvsForLeafStorage(vc);
			break;
		}
	case COMPRESSED:
		{
			pvsSize = GetEntriesInPvsForCompressedStorage(vc);
			break;
		}
	case PVS_IN_INTERIORS:
	default:
		// pvs is stored consistently in the tree up to the root
		// just return pvs size 
		pvsSize = vc->GetPvs().GetSize();	
		break;
	}

	return pvsSize;  
}


float ViewCellsTree::GetMemoryCost(ViewCell *vc) const
{
	return (float)CountStoredPvsEntries(vc) * ObjectPvs::GetEntrySize();
}

//#if HAS_TO_BE_REDONE
int ViewCellsTree::CountStoredPvsEntries(ViewCell *root) const
{
	int pvsSize = root->GetPvs().GetSize();

	// recursivly count leaves
	if (!root->IsLeaf())
	{
		ViewCellInterior *interior = static_cast<ViewCellInterior *>(root);

		ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

		for (it = interior->mChildren.begin(); it != it_end; ++ it)
		{
			pvsSize += CountStoredPvsEntries(*it);
		}
	}

	return pvsSize;    	
}
//#endif

int ViewCellsTree::ViewCellsStorage() const
{
	return mViewCellsStorage;
}


ViewCell *ViewCellsTree::GetActiveViewCell(ViewCellLeaf *vc) const
{
	return vc->GetActiveViewCell();
}


void ViewCellsTree::PropagatePvs(ViewCell *root)
{	
	ViewCell *viewCell = root;

	// propagate pvs up
	while (viewCell->GetParent())
	{
		ObjectPvs mergedPvs;
		viewCell->GetParent()->GetPvs().MergeInPlace(root->GetPvs());

		viewCell = viewCell->GetParent();
	}

	if (root->IsLeaf())
		return;

	// propagate pvs to the leaves
	stack<ViewCell *> tstack;
	tstack.push(root);

	while (!tstack.empty())
	{
		ViewCell *viewCell = tstack.top();
		tstack.pop();

		if (viewCell != root)
		{
			viewCell->GetPvs().MergeInPlace(root->GetPvs());
		}

		if (!viewCell->IsLeaf())
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(viewCell);

			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}
		}
	}
}


void ViewCellsTree::AssignRandomColors()
{
	TraversalQueue tqueue;
	tqueue.push(mRoot);
	
	mRoot->SetColor(RandomColor(0.3f, 1.0f));
	
	while (!tqueue.empty())
	{
		ViewCell *vc = tqueue.top();
		tqueue.pop();

		// save the view cells if it is a leaf or if enough view cells 
		// have already been traversed because of the priority queue, 
		// this will be the optimal set of v
		if (!vc->IsLeaf()) 
		{	
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
		  
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
		  
			float maxProbability = -1.0f;
		  
			ViewCell *maxViewCell = NULL;
			for (it = interior->mChildren.begin(); it != it_end; ++ it) 
			{
				ViewCell *v = *it;
			  
				// set random color
				v->SetColor(RandomColor(0.3f, 1.0f));
			  
				if (v->GetVolume() > maxProbability) 
				{
					maxProbability = v->GetVolume();
					maxViewCell = v;
				}

				if (maxViewCell)
				{
					maxViewCell->SetColor(vc->GetColor());
				}
				
				tqueue.push(v);
			}
		}	 
	}
}


/** Get costs resulting from each merge step. 
*/
void ViewCellsTree::GetCostFunction(vector<float> &costFunction)
{
	TraversalQueue tqueue;
	tqueue.push(mRoot);

	int numViewCells = 1;

	const float vol = mViewCellsManager->GetViewSpaceBox().GetVolume();
	const float rootPvs = GetPvsCost(mRoot);
	
	float totalRenderCost;

	float totalPvsCost = rootPvs;
	totalRenderCost = (float)rootPvs;

	costFunction.push_back(totalRenderCost);

	//-- go through tree in the order of render cost decrease
	//-- which is the same order as the view cells were merged
	//-- or the reverse order of subdivision for subdivision-only 
	//-- view cell hierarchies.

	while (!tqueue.empty())
	{
		ViewCell *vc = tqueue.top();
		tqueue.pop();

		if (!vc->IsLeaf()) 
		{	
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

			const float parentCost = GetPvsCost(interior);
			const float parentExpCost = parentCost * interior->GetVolume();

			-- numViewCells;

			float childExpCost = 0;
			float childCost = 0;
			
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				ViewCell *vc = *it;

				const float pvsCost = GetPvsCost(vc);
				
				childExpCost += (float) pvsCost * vc->GetVolume();
				childCost += pvsCost;
				
				tqueue.push(vc);
				++ numViewCells;
			}

			// update stats for this view cell
			const float costDecr = (parentExpCost - childExpCost) / vol;
			totalRenderCost -= costDecr;
			
			costFunction.push_back(totalRenderCost);
		}
	}
}


/** Get storage costs resulting from each merge step. 
*/
void ViewCellsTree::GetStorageFunction(vector<int> &storageFunction)
{
	TraversalQueue tqueue;
	tqueue.push(mRoot);

	int numViewCells = 1;

	const float vol = mViewCellsManager->GetViewSpaceBox().GetVolume();
	const int rootEntries = GetPvsEntries(mRoot);

	int entriesInPvs = rootEntries;
	// one entry into the pvs
    const int entryStorage = sizeof(PvsData) + sizeof(int); 

	storageFunction.push_back(rootEntries);

	////////////
	//-- go through tree in the order of render cost decrease
	//-- which is the same order as the view cells were merged
	//-- or the reverse order of subdivision for subdivision-only 
	//-- view cell hierarchies.

	while (!tqueue.empty())
	{
		ViewCell *vc = tqueue.top();
		tqueue.pop();

		if (!vc->IsLeaf()) 
		{	
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

			const float parentPvsCost = GetPvsCost(interior);
			const int parentPvsEntries = GetPvsEntries(interior);
            const float parentExpCost = (float)parentPvsCost * interior->GetVolume();

			float childExpCost = 0;
			int childPvs = 0;
			int childPvsEntries = 0;

			-- numViewCells;

			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				ViewCell *vc = *it;

				const int pvsEntries = GetPvsEntries(vc);
				childPvsEntries += pvsEntries;

				tqueue.push(vc);
				++ numViewCells;
			}

			// update stats for this view cell
			const float costDecr = (parentExpCost - childExpCost) / vol;

			entriesInPvs += childPvsEntries - parentPvsEntries;

			const int storageCost = entriesInPvs * entryStorage;
			storageFunction.push_back(storageCost);
		}
	}
}


void ViewCellsTree::UpdateViewCellsStats(ViewCell *vc, 
										 ViewCellsStatistics &vcStat)
{
	++ vcStat.viewCells;
		
	const float pvsCost = GetPvsCost(vc);

	vcStat.pvsCost += pvsCost;

	if (pvsCost < Limits::Small)
		++ vcStat.emptyPvs;

	if (pvsCost > vcStat.maxPvs)
		vcStat.maxPvs = pvsCost;

	if (pvsCost < vcStat.minPvs)
		vcStat.minPvs = pvsCost;

	if (!vc->GetValid())
		++ vcStat.invalid;
}


bool ViewCellsTree::Export(OUT_STREAM &stream, const bool exportPvs)
{
	// export recursivly all view cells from the root
	ExportViewCell(mRoot, stream, exportPvs);

	return true;
}


void ViewCellsTree::CreateUniqueViewCellsIds()
{
	stack<ViewCell *> tstack;

	int currentId = 0;

	tstack.push(mRoot);

	while (!tstack.empty())
	{
		ViewCell *vc = tstack.top();
		tstack.pop();

		if (vc->GetId() != -1) // out of bounds
			vc->SetId(currentId ++);

		if (!vc->IsLeaf())
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
			
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}
		}
	}
}


void ViewCellsTree::ExportPvs(ViewCell *viewCell, OUT_STREAM &stream)
{
	ObjectPvsIterator it = viewCell->GetPvs().GetIterator();

	while (it.HasMoreEntries())
	{
		Intersectable *obj = it.Next();

		// hack: just output all the "elementary" objects
		if (0 && (obj->Type() == Intersectable::BVH_INTERSECTABLE))
		{
			ObjectContainer objects; 
			BvhNode *node = static_cast<BvhNode *>(obj);
			node->CollectObjects(objects);
		
			ObjectContainer::const_iterator oit, oit_end = objects.end();
			for (oit = objects.begin(); oit != oit_end; ++ oit)
			{
				stream << (*oit)->GetId() << " ";
			}
		}
		else
		{
			stream << obj->GetId() << " ";
		}
	}
}


void ViewCellsTree::ExportViewCell(ViewCell *viewCell, 
								   OUT_STREAM &stream, 
								   const bool exportPvs)
{
	if (viewCell->IsLeaf())
	{
		stream << "<Leaf ";
		stream << "id=\"" << viewCell->GetId() << "\" ";
		stream << "active=\"" << static_cast<ViewCellLeaf *>(viewCell)->GetActiveViewCell()->GetId() << "\" ";
		stream << "mergecost=\"" << viewCell->GetMergeCost() << "\" ";
		stream << "pvs=\"";
		
		//-- export pvs, i.e., the ids of the objects in the pvs
		if (exportPvs)
		{
			ExportPvs(viewCell, stream);
		}
		stream << "\" />" << endl;
	}
	else
	{
		ViewCellInterior *interior = static_cast<ViewCellInterior *>(viewCell);
	
		stream << "<Interior ";
		stream << "id=\"" << viewCell->GetId() << "\" ";
		stream << "mergecost=\"" << viewCell->GetMergeCost() << "\" ";
		stream << "pvs=\"";

		// NOTE: do not export pvss for interior view cells because
		// they can be completely reconstructed from the leaf pvss
		// on the other hand: we could store a tag with the compression scheme,
        // then some scheme were pvs is in the interiors could be used
		if (exportPvs)
		{
			ExportPvs(viewCell, stream);
		}

		stream << "\" >" << endl;

		//-- recursivly export child view cells
		ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

		for (it = interior->mChildren.begin(); it != it_end; ++ it)
		{
			ExportViewCell(*it, stream, exportPvs);
		}

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


void ViewCellsTree::ResetPvs()
{
	stack<ViewCell *> tstack;

	tstack.push(mRoot);

	while (!tstack.empty())
	{
		ViewCell *vc = tstack.top();
		tstack.pop();

		vc->GetPvs().Clear();
		
		if (!vc->IsLeaf())
		{
			ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
			ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

			for (it = interior->mChildren.begin(); it != it_end; ++ it)
			{
				tstack.push(*it);
			}
		}
	}
}


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


void ViewCellsTree::SetActiveSetToLeaves()
{
	// todo
}



/**************************************************************************/
/*                     MergeCandidate implementation                      */
/**************************************************************************/


MergeCandidate::MergeCandidate(ViewCell *l, ViewCell *r):
mRenderCost(0),
mDeviationIncr(0),
mLeftViewCell(l),
mRightViewCell(r),
mInitialLeftViewCell(l),
mInitialRightViewCell(r)
{
	//EvalMergeCost();
}


void MergeCandidate::SetRightViewCell(ViewCell *v)
{
	mRightViewCell = v;
}


void MergeCandidate::SetLeftViewCell(ViewCell *v)
{
	mLeftViewCell = v;
}


ViewCell *MergeCandidate::GetRightViewCell() const
{
	return mRightViewCell;
}


ViewCell *MergeCandidate::GetLeftViewCell() const
{
	return mLeftViewCell;
}


ViewCell *MergeCandidate::GetInitialRightViewCell() const
{
	return mInitialRightViewCell;
}


ViewCell *MergeCandidate::GetInitialLeftViewCell() const
{
	return mInitialLeftViewCell;
}


bool MergeCandidate::IsValid() const
{
	// if one has a parent, it was already merged
	return !(mLeftViewCell->GetParent() || mRightViewCell->GetParent());
}


float MergeCandidate::GetRenderCost() const
{
	return mRenderCost;
}


float MergeCandidate::GetDeviationIncr() const
{
	return mDeviationIncr;
}


float MergeCandidate::GetMergeCost() const
{
	return mRenderCost * sRenderCostWeight + 
		   mDeviationIncr * (1.0f - sRenderCostWeight);
}



/************************************************************************/
/*                    MergeStatistics implementation                    */
/************************************************************************/


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

	app << setprecision(4);

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

	app << "#N_CCTIME ( Collect candidates time [s] )\n" << collectTime * 1e-3f << " \n";

	app << "#N_MTIME ( Merge time [s] )\n" << mergeTime * 1e-3f << " \n";

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

	app << "#N_CANDIDATES ( Number of merge candidates )\n" << candidates << "\n";

	app << "#N_MERGEDSIBLINGS ( Number of merged siblings )\n" << siblings << "\n";

	app << "#OVERALLCOST ( overall merge cost )\n" << overallCost << "\n";

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

	app << "#MAX_TREEDIST ( Maximal distance in tree of merged leaves )\n" << maxTreeDist << "\n";

	app << "#AVG_TREEDIST ( Average distance in tree of merged leaves )\n" << AvgTreeDist() << "\n";

	app << "#EXPECTEDCOST ( expected render cost )\n" << expectedRenderCost << "\n";

	app << "#DEVIATION ( deviation )\n" << deviation << "\n";

	app << "#HEURISTICS ( heuristics )\n" << heuristics << "\n";
	

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


}