// ================================================================
// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
// ****************************************************************
/**
   The KD tree based LSDS
*/
// Initial coding by
/**
   @author Jiri Bittner
*/

// Standard headers
#include <stack>
#include <queue>
#include <algorithm>
#include <fstream>
#include <string>

#include "RssTree.h"

#include "Environment.h"
#include "VssRay.h"
#include "Intersectable.h"
#include "Ray.h"
#include "Containers.h"
#include "ViewCell.h"
#include "Exporter.h"
#include "Preprocessor.h"
#include "SceneGraph.h"
#include "SamplingStrategy.h"
#include "ViewCellsManager.h"


namespace GtpVisibilityPreprocessor {


#define DEBUG_SPLIT_COST 0
#define DEBUG_SPLITS 0

#define EVAL_VIEWCELLS 0
  
// Static variables
int
RssTreeLeaf::mailID = 0;

inline void
AddObject2Pvs(Intersectable *object,
			  const int side,
			  int &pvsBack,
			  int &pvsFront)
{
  
  if (!object)
	return;
  
  if (side <= 0) {
	if (!object->Mailed() && !object->Mailed(2)) {
	  pvsBack++;
	  if (object->Mailed(1))
		object->Mail(2);
	  else
		object->Mail();
	}
  }
  
  if (side >= 0) {
	if (!object->Mailed(1) && !object->Mailed(2)) {
	  pvsFront++;
	  if (object->Mailed())
		object->Mail(2);
	  else
		object->Mail(1);
	}
  }
}

inline void
AddViewcells2Pvs(const ViewCellContainer &viewcells,
				 const int side,
				 int &viewcellsBack,
				 int &viewcellsFront)
{
  ViewCellContainer::const_iterator it = viewcells.begin();
  
  for (; it != viewcells.end(); ++it) {
	ViewCell *viewcell = *it;
	if (side <= 0) {
	  if (!viewcell->Mailed() && !viewcell->Mailed(2)) {
		viewcellsBack++;
		if (viewcell->Mailed(1))
		  viewcell->Mail(2);
		else
		  viewcell->Mail();
	  }
	}
	
	if (side >= 0) {
	  if (!viewcell->Mailed(1) && !viewcell->Mailed(2)) {
		viewcellsFront++;
		if (viewcell->Mailed())
		  viewcell->Mail(2);
		else
		  viewcell->Mail(1);
	  }
	}
  }
}


// Constructor
RssTree::RssTree()
{
  Environment::GetSingleton()->GetIntValue("RssTree.maxDepth", termMaxDepth);
  Environment::GetSingleton()->GetIntValue("RssTree.minPvs", termMinPvs);
  Environment::GetSingleton()->GetIntValue("RssTree.minRays", termMinRays);
  Environment::GetSingleton()->GetFloatValue("RssTree.maxRayContribution", termMaxRayContribution);
  Environment::GetSingleton()->GetFloatValue("RssTree.maxCostRatio", termMaxCostRatio);

  Environment::GetSingleton()->GetFloatValue("RssTree.minSize", termMinSize);
  termMinSize = sqr(termMinSize);
	
  Environment::GetSingleton()->GetFloatValue("RssTree.refDirBoxMaxSize", refDirBoxMaxSize);
  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
  
  Environment::GetSingleton()->GetFloatValue("RssTree.epsilon", epsilon);
  Environment::GetSingleton()->GetFloatValue("RssTree.ct_div_ci", ct_div_ci);
	
  Environment::GetSingleton()->GetFloatValue("RssTree.maxTotalMemory", maxTotalMemory);
  Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
  
  Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);


  
  Environment::GetSingleton()->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
  //= 1000;
  Environment::GetSingleton()->GetIntValue("RssTree.minCollapseDepth", minCollapseDepth);
  //  int minCollapseDepth = 4;

  //  pRefDirThresh = cos(0.5*M_PI - M_PI*refDirAngle/180.0);
  //  cosRefDir = cos(M_PI*refDirAngle/180.0);
  //  sinRefDir = sin(M_PI*refDirAngle/180.0);
  
  
  // split type
  char sname[128];
  Environment::GetSingleton()->GetStringValue("RssTree.splitType", sname);
  string name(sname);
	
  if (name.compare("regular") == 0)
    splitType = ESplitRegular;
  else
    if (name.compare("heuristic") == 0)
      splitType = ESplitHeuristic;
	else
	  if (name.compare("hybrid") == 0)
		splitType = ESplitHybrid;
	  else {
		cerr<<"Invalid RssTree split type "<<name<<endl;
		exit(1);
	  }

  Environment::GetSingleton()->GetIntValue("RssTree.hybridDepth", mHybridDepth);

  Environment::GetSingleton()->GetBoolValue("RssTree.randomize", randomize);
  Environment::GetSingleton()->GetBoolValue("RssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);

  Environment::GetSingleton()->GetBoolValue("RssTree.interleaveDirSplits", mInterleaveDirSplits);
  Environment::GetSingleton()->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);

  Environment::GetSingleton()->GetBoolValue("RssTree.importanceBasedCost", mImportanceBasedCost);
  
  Environment::GetSingleton()->GetIntValue("RssTree.maxRays", mMaxRays);

  Environment::GetSingleton()->GetBoolValue("RssTree.perObjectTree", mPerObjectTree);
  
  //  mRoots;
  
  splitCandidates = new vector<SortableEntry>;
}


RssTree::~RssTree()
{
  for (int i=0; i < mRoots.size(); i++)
	if (mRoots[i])
	  delete mRoots[i];
}




void
RssStatistics::Print(ostream &app) const
{
  app << "###### RssTree statistics ######\n";

  app << "#N_RAYS ( Number of rays )\n"
      << rays <<endl;

  app << "#N_INITPVS ( Initial PVS size )\n"
      << initialPvsSize <<endl;
  
  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";

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

  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
  for (int i=0; i<7; i++)
    app << splits[i] <<" ";
  app <<endl;

  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
    rayRefs << "\n";

  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
    rayRefs/(double)rays << "\n";

  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
    rayRefs/(double)Leaves() << "\n";

  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
    maxRayRefs << "\n";


  //  app << setprecision(4);

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

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

  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minRays )\n"<<
    minRaysNodes*100/(double)Leaves()<<endl;
	
  app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
    minSizeNodes*100/(double)Leaves()<<endl;

  app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
    maxRayContribNodes*100/(double)Leaves()<<endl;

  app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
    maxCostRatioNodes*100/(double)Leaves()<<endl;

  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
    addedRayRefs<<endl;

  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
    removedRayRefs<<endl;

  //  app << setprecision(4);

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

  app << "###### END OF RssTree statistics ######\n";

}


void
RssTree::UpdatePvsSize(RssTreeLeaf *leaf)
{
  if (!leaf->mValidPvs) {
	Intersectable::NewMail();
	int pvsSize = 0;
	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		ri != leaf->rays.end();
		ri++)
	  if ((*ri).mRay->IsActive()) {
		Intersectable *object;
		object = (*ri).GetObject();
		if (object && !object->Mailed()) {
		  pvsSize++;
		  object->Mail();
		}
	  }
	leaf->SetPvsSize(pvsSize);
	
	ComputeImportance(leaf);
  }
}

bool
RssTree::ClipRay(
				 RssTreeNode::RayInfo &rayInfo,
				 const AxisAlignedBox3 &box
				 )
{
  float tmin, tmax;
  static Ray ray;
  cerr<<"Clip not reimplented yet...Quiting\n";
  exit(1);
  ray.Init(rayInfo.GetOrigin(), rayInfo.GetDir(), Ray::LINE_SEGMENT);

  if (!box.ComputeMinMaxT(ray, &tmin, &tmax) || tmin>=tmax)
	return false;
  
  // now check if the ray origin lies inside the box
  if ( tmax < rayInfo.mRay->GetSize() ) {
	// this ray does not leave the box
	rayInfo.mRay->SetupEndPoints(
								 ray.Extrap(tmax),
								 rayInfo.mRay->mTermination
								 );
	return true;
  }

  return false;
}


void
RssTree::Construct(
				   ObjectContainer &objects,
				   VssRayContainer &rays
				   // forced bounding box is only used when computing from-box
				   // visibility
				   //				   AxisAlignedBox3 *forcedBoundingBox
				   )
{
  cout<<"Constructing rss tree"<<endl<<flush;
  stat.Start();
  
  maxMemory = maxStaticMemory;

  //  if (root)
  //    delete root;

  if (mPerObjectTree) {
	// get max id from the rays
	int i;
	mRoots.resize(objects.size());
	for (i = 0; i < objects.size(); i++) {
	  RssTreeLeaf *leaf = new RssTreeLeaf(NULL, 0);
	  //	  leaf->bbox.Initialize();
	  leaf->dirBBox.Initialize();
	  leaf->dirBBox.SetMin(2, 0.0f);
	  leaf->dirBBox.SetMax(2, 0.0f);
	  mRoots[i] = leaf;
	}

	// init spatial bounding boxes
	for (i = 0; i < objects.size(); i++) {
	  GetRoot(objects[i])->bbox = objects[i]->GetBox();
	}
	stat.nodes = i;
	stat.leaves = i;
  } else {
	mRoots.resize(1);
	RssTreeLeaf *leaf = new RssTreeLeaf(NULL, (int)rays.size());
	leaf->bbox.Initialize();
	leaf->dirBBox.Initialize();
	leaf->dirBBox.SetMin(2, 0.0f);
	leaf->dirBBox.SetMax(2, 0.0f);
	mRoots[0] = leaf;
	stat.nodes = 1;
	stat.leaves = 1;
  }

  for(VssRayContainer::const_iterator ri = rays.begin();
	  ri != rays.end();
	  ri++) {
	
	  RssTreeNode::RayInfo info(*ri);
	  if (!EpsilonEqualV3(info.GetOrigin(), info.GetTermination(), Limits::Small)) {

	  // first construct a leaf that will get subdivide
	  RssTreeLeaf *leaf = (RssTreeLeaf *) GetRoot(info.GetSourceObject());
	  

	  leaf->AddRay(info);
	  
	  // leaf bbox contains bbox of origins only

	  // include both origin and terminatin in the global bbox
#if USE_ORIGIN
	  leaf->bbox.Include((*ri)->GetOrigin());
	  bbox.Include((*ri)->GetOrigin());
#else
	  bbox.Include((*ri)->GetTermination());
	  leaf->bbox.Include((*ri)->GetTermination());
#endif	  
	  Vector3 dVec = Vector3(
							 (*ri)->GetDirParametrization(0),
							 (*ri)->GetDirParametrization(1),
							 0
							 );

	  leaf->dirBBox.Include(dVec);
	  dirBBox.Include(dVec);
	  }
  }

  // make the z axis (unused) a unit size
  // important for volume computation
  
  //  if ( forcedBoundingBox ) 
  //	bbox = *forcedBoundingBox;
  
  cout<<"Bbox = "<<bbox<<endl;
  cout<<"Dirr Bbox = "<<dirBBox<<endl;
  
  stat.rays = rays.size();
  stat.initialPvsSize = 0;
  for (int i=0; i < mRoots.size(); i++) {
	RssTreeLeaf *leaf = (RssTreeLeaf *)mRoots[i];
	UpdatePvsSize(leaf);
	stat.initialPvsSize += leaf->GetPvsSize();
	mRoots[i] = Subdivide(TraversalData(leaf, GetBBox(leaf), 0));
  }
  
  if (splitCandidates) {
    // force realease of this vector
    delete splitCandidates;
    splitCandidates = new vector<SortableEntry>;
  }
  
  stat.Stop();
  stat.Print(cout);
  cout<<"#Total memory="<<GetMemUsage()<<endl;

  // this rotine also updates importances etc...
}

int
RssTree::UpdateSubdivision()
{
  priority_queue<TraversalData> tStack;
  //  stack<TraversalData> tStack;
  
  //  tStack.push(TraversalData(root, bbox, 0));
  PushRoots(tStack);
			
  AxisAlignedBox3 backBox;
  AxisAlignedBox3 frontBox;

  maxMemory = maxTotalMemory;
  int subdivided = 0;
  int lastMem = 0;
  while (!tStack.empty()) {
		
	float mem = GetMemUsage();
		
	if ( lastMem/10 != ((int)mem)/10) {
	  cout<<mem<<" MB"<<endl;
	}
	lastMem = (int)mem;
		
	if (  mem > maxMemory ) {
      // count statistics on unprocessed leafs
      while (!tStack.empty()) {
		//				EvaluateLeafStats(tStack.top());
		tStack.pop();
      }
      break;
    }
    
    TraversalData data = tStack.top();
    tStack.pop();

	if (data.node->IsLeaf()) {
	  RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
										data.bbox,
										backBox,
										frontBox
										);
	  if (!node->IsLeaf()) {
		subdivided++;
		
		
		RssTreeInterior *interior = (RssTreeInterior *) node;
		// push the children on the stack
		tStack.push(TraversalData(interior->back, backBox, data.depth+1));
		tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
	  } else {
		//      EvaluateLeafStats(data);
	  }
	} else {
	  RssTreeInterior *interior = (RssTreeInterior *) data.node;
	  tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
	  tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
	}
  }
  return subdivided;
}


RssTreeNode *
RssTree::Subdivide(const TraversalData &tdata)
{
  RssTreeNode *result = NULL;

  priority_queue<TraversalData> tStack;
  //  stack<TraversalData> tStack;
  
  tStack.push(tdata);

  AxisAlignedBox3 backBox;
  AxisAlignedBox3 frontBox;

  
  int lastMem = 0;
  while (!tStack.empty()) {

	float mem = GetMemUsage();
		
	if ( lastMem/10 != ((int)mem)/10) {
	  cout<<mem<<" MB"<<endl;
	}
	lastMem = (int)mem;
		
	if (  mem > maxMemory ) {
      // count statistics on unprocessed leafs
      while (!tStack.empty()) {
		EvaluateLeafStats(tStack.top());
		tStack.pop();
      }
      break;
    }
    
    TraversalData data = tStack.top();
    tStack.pop();

#if DEBUG_SPLITS
	Debug<<"#Splitting node"<<endl;
	data.node->Print(Debug);
#endif
	RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
									  data.bbox,
									  backBox,
									  frontBox
									  );
	

	if (result == NULL)
      result = node;
    
    if (!node->IsLeaf()) {
			
      RssTreeInterior *interior = (RssTreeInterior *) node;
      // push the children on the stack
      tStack.push(TraversalData(interior->back, backBox, data.depth+1));
      tStack.push(TraversalData(interior->front, frontBox, data.depth+1));

#if DEBUG_SPLITS
	  Debug<<"#New nodes"<<endl;
	  interior->back->Print(Debug);
	  interior->front->Print(Debug);
	  Debug<<"#####################################"<<endl;
#endif
	  
    } else {
      EvaluateLeafStats(data);
    }
  }

  return result;
}


// returns selected plane for subdivision
int
RssTree::SelectPlane(
					 RssTreeLeaf *leaf,
					 const AxisAlignedBox3 &box,
					 SplitInfo &info
					 )
{

  BestCostRatio(leaf,
				info);
  
#if DEBUG_SPLIT_COST
  Debug<<"Split Info:"<<endl;
  Debug<<"axis="<<info.axis<<" ratio="<<info.costRatio<<endl;
  Debug<<"rays="<<info.rays<<
  	" rays back="<<info.raysBack<<
  	" rays front="<<info.raysFront<<endl;
  Debug<<"c="<<info.contribution<<
  	" c back="<<info.contributionBack<<
  	" c front="<<info.contributionFront<<endl;

  //  Debug<<"viewcells="<<info.viewCells<<
  //	" viewcells back="<<info.viewCellsBack<<
  //	" viewcells back="<<info.viewCellsFront<<endl;
#endif
  
  if (info.costRatio > termMaxCostRatio) {
	//		cout<<"Too big cost ratio "<<costRatio<<endl;
	stat.maxCostRatioNodes++;
	return -1;
  }
	
#if 0
  cout<<
	"pvs="<<leaf->mPvsSize<<
	" rays="<<leaf->rays.size()<<
	" rc="<<leaf->GetAvgRayContribution()<<
	" axis="<<info.axis<<endl;
#endif
  
  return info.axis;
}


void
RssTree::GetCostRatio(
					  RssTreeLeaf *leaf,
					  SplitInfo &info
					  )
{

  AxisAlignedBox3 box;
  float minBox, maxBox;

  if (info.axis < 3) {
	box = GetBBox(leaf);
	minBox = box.Min(info.axis);
	maxBox = box.Max(info.axis);
  }	else {
	box = GetDirBBox(leaf);
	minBox = box.Min(info.axis-3);
	maxBox = box.Max(info.axis-3);
  }
	
  float sizeBox = maxBox - minBox;

  int pvsSize = leaf->GetPvsSize();

  if (!mImportanceBasedCost) {
	const int costMethod = 0;
	
	switch (costMethod) {
	case 0: {
	  float sum = info.pvsBack*(info.position - minBox) + info.pvsFront*(maxBox - info.position);
	  float newCost = ct_div_ci + sum/sizeBox;
	  float oldCost = (float)pvsSize;
	  info.costRatio = newCost/oldCost;
	  break;
	}
	case 6: {
	  float sum = info.raysBack*(info.position - minBox) + info.raysFront*(maxBox - info.position);
	  float newCost = ct_div_ci + sum/sizeBox;
	  float oldCost = (float)info.rays;
	  info.costRatio = newCost/oldCost;
	  break;
	}
	case 7: {
	  float benefit =
		(sqr(info.contributionBack) + sqr(info.contributionFront))/
		(2.0f*sqr(info.contribution));
	  info.costRatio = 1.0f/(benefit + Limits::Small);
	  break;
	}

	  
	case 1: {
	  float newContrib =
		(info.contributionBack*(info.position - minBox) + 
		 +
		 info.contributionFront*(maxBox - info.position))/sizeBox;
	  float oldContrib = info.contribution;
	  info.costRatio = newContrib/oldContrib;
	  //	  cout<<info.contribution<<" "<<info.contributionBack<<" "<<info.contributionFront<<endl;
	  break;
	} 
	case 2: {
	  float sum =
		info.viewCellsBack*(info.position - minBox) +
		info.viewCellsFront*(maxBox - info.position);
	  float newCost = ct_div_ci + sum/sizeBox;
	  float oldCost = (float)info.viewCells;
	  info.costRatio = newCost/oldCost;
	  break;
	}
	case 3: {
	  float newCost = (float)(info.raysBack*info.pvsBack  + info.raysFront*info.pvsFront);
	  float oldCost = (float)leaf->rays.size()*pvsSize;
	  info.costRatio = newCost/oldCost;
	  break;
	}
	}
  } else {
	const int costMethod = 1;
	// importance based cost
	switch (costMethod) {
	case 0: {
	  break;
	}
	case 1: {
	  float newContrib =
		sqr(info.pvsBack/(info.raysBack + Limits::Small)) +
		sqr(info.pvsFront/(info.raysFront + Limits::Small));
	  float oldContrib = sqr(leaf->GetAvgRayContribution());
	  info.costRatio = oldContrib/newContrib;
	  break;
	}
	case 2: {
	  float newCost = (info.pvsBack  + info.pvsFront)*0.5f;
	  float oldCost = (float)pvsSize;
	  info.costRatio = newCost/oldCost;
	  break;
	}
	case 3: {
	  float newCost = (float)abs(info.raysBack  - info.raysFront);
	  float oldCost = (float)leaf->rays.size();
	  info.costRatio = newCost/oldCost;
	  break;
	}
  }
  }
}
							

void
RssTree::EvalCostRatio(
					   RssTreeLeaf *leaf,
					   SplitInfo &info
					   )
{
  info.rays = 0;
  info.raysBack = 0;
  info.raysFront = 0;
  info.pvsFront = 0;
  info.pvsBack = 0;
  info.viewCells = 0;
  info.viewCellsFront = 0;
  info.viewCellsBack = 0;



  float sumWeights = Limits::Small;
  float sumWeightsBack = Limits::Small;
  float sumWeightsFront = Limits::Small;
  
  float sumContribution = 0.0f;
  float sumContributionBack = 0.0f;
  float sumContributionFront = 0.0f;

  Intersectable::NewMail();

#if EVAL_VIEWCELLS
  // count the numebr of viewcells first
  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	  ri != leaf->rays.end();
	  ri++)
	if ((*ri).mRay->IsActive()) {

	  ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
	  for (; it != (*ri).mRay->mViewCells.end(); ++it) {
		if (!(*it)->Mailed()) {
		  (*it)->Mail();
		  info.viewCells++;
		}
	  }
	}
#endif
  
  Intersectable::NewMail(3);
  
  // this is the main ray classification loop!
  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	  ri != leaf->rays.end();
	  ri++)
	if ((*ri).mRay->IsActive()) {
	  int side;
	  
	  
	  // determine the side of this ray with respect to the plane
	  side = (*ri).ComputeRaySide(info.axis, info.position);
	  
	  float weight, contribution;
	  
	  GetRayContribution(*ri, weight, contribution);
	  sumWeights += weight;
	  sumContribution += contribution;
	  info.rays++;
	  if (side <= 0) {
		info.raysBack++;
		sumWeightsBack += weight;
		sumContributionBack += contribution;
	  }
	  
	  if (side >= 0) {
		info.raysFront++;
		sumWeightsFront += weight;
		sumContributionFront += contribution;
	  }
	  
	  AddObject2Pvs((*ri).GetObject(), side, info.pvsBack, info.pvsFront);

#if EVAL_VIEWCELLS
	  AddViewcells2Pvs((*ri).mRay->mViewCells,
					   side,
					   info.viewCellsBack,
					   info.viewCellsFront);
#endif
	}

  if (sumWeights!=0.0f)
	info.contribution = sumContribution/sumWeights;
  else
	info.contribution = 0.0f;

  if (sumWeightsBack!=0.0f)
	info.contributionBack = sumContributionBack/sumWeightsBack;
  else
	info.contributionBack = 0.0f;

  if (sumWeightsFront!=0.0f)
	info.contributionFront = sumContributionFront/sumWeightsFront;
  else
	info.contributionFront = 0.0f;

  //  info.costRatio = 0.1 + Random(0.5f);
  //  return;

  GetCostRatio(
			   leaf,
			   info);
  
  //cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
  //  float oldCost = leaf->rays.size();
  
  //  cout<<"ratio="<<ratio<<endl;
}

void
RssTree::BestCostRatio(
					   RssTreeLeaf *leaf,
					   SplitInfo &info
					   )
{
  SplitInfo nInfo[5];
  int bestAxis = -1;
  
  AxisAlignedBox3 sBox = GetBBox(leaf);
  AxisAlignedBox3 dBox = GetDirBBox(leaf);
  // int sAxis = box.Size().DrivingAxis();
  int sAxis = sBox.Size().DrivingAxis();
  int dAxis = dBox.Size().DrivingAxis() + 3;
  
  float dirSplitBoxSize = 0.01f;
  bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
		
  bool useCyclingAxis = true;

  if (splitType == ESplitHeuristic ||
	  (splitType == ESplitHybrid && leaf->depth > mHybridDepth))
	useCyclingAxis = false;
  
  int cyclingAxis = 0;
  
  if (leaf->parent) 
	cyclingAxis = (leaf->parent->axis + 1)%5;
  
  if (0 && useCyclingAxis) {
	int axis = cyclingAxis;
	info.position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
	info.axis = axis;
	info.costRatio = 0.5f; // not true but good for speeding up the subdivision at the top levels!
	info.raysBack = info.raysFront = leaf->rays.size()/2;
	return;
  }

  
  for (int axis = 0; axis < 5; axis++)
	if (
		(axis < 3 && (leaf->depth < mDirSplitDepth ||  mInterleaveDirSplits)) ||
		(axis >= 3 && (leaf->depth >= mDirSplitDepth))
		) {
	  nInfo[axis].axis = axis;
	  if ((!useCyclingAxis && (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis))
		  || (useCyclingAxis && (axis == cyclingAxis))
		  ) {
		
		if (splitType == ESplitRegular) {
		  if (axis < 3)
			nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
		  else
			nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
		  EvalCostRatio(leaf,
						nInfo[axis]);
		} else
		  if (splitType == ESplitHeuristic) {
			EvalCostRatioHeuristic(
								   leaf,
								   nInfo[axis]
								   );
		  } else
			if (splitType == ESplitHybrid) {
			  if (leaf->depth > mHybridDepth)
				EvalCostRatioHeuristic(
									   leaf,
									   nInfo[axis]
									   );
			  else {
				if (axis < 3)
				  nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
				else
				  nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
				
				EvalCostRatio(leaf,
							  nInfo[axis]
							  );
			  }
			} else {
			  cerr<<"RssTree: Unknown split heuristics\n";
			  exit(1);
			}
		
		if ( bestAxis == -1)
		  bestAxis = axis;
		else
		  if ( nInfo[axis].costRatio < nInfo[bestAxis].costRatio )
			bestAxis = axis;
	  }
	}
  
  info = nInfo[bestAxis];
}

	
void
RssTree::EvalCostRatioHeuristic(
								RssTreeLeaf *leaf,
								SplitInfo &info
								)
{
  AxisAlignedBox3 box;
  float minBox, maxBox;
	
  if (info.axis < 3) {
	box = GetBBox(leaf);
	minBox = box.Min(info.axis);
	maxBox = box.Max(info.axis);
  } else {
	box = GetDirBBox(leaf);
	minBox = box.Min(info.axis-3);
	maxBox = box.Max(info.axis-3);
  }
	
  SortSubdivisionCandidates(leaf, info.axis);
  
  // go through the lists, count the number of objects left and right
  // and evaluate the following cost funcion:
  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
  
  SplitInfo currInfo;
  currInfo.axis = info.axis;

  currInfo.raysBack = 0;
  currInfo.raysFront = leaf->rays.size();
  
  currInfo.pvsBack = 0;
  currInfo.pvsFront = leaf->GetPvsSize();


  float sumWeights = Limits::Small;
  float sumWeightsBack = Limits::Small;
  float sumWeightsFront = Limits::Small;
  
  float sumContribution = 0.0f;
  float sumContributionBack = 0.0f;
  float sumContributionFront = 0.0f;

  
  float sizeBox = maxBox - minBox;
  
  float minBand = minBox + 0.1f*(maxBox - minBox);
  float maxBand = minBox + 0.9f*(maxBox - minBox);
	
  // best cost ratio
  info.costRatio = 1e20f;

  currInfo.viewCells = 0;
  currInfo.rays = 0;

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

	  float weight, contribution;
	  GetRayContribution(*ri, weight, contribution);

	  sumWeights += weight;
	  sumContribution += contribution;
	  currInfo.rays++;

	  Intersectable *object = (*ri).GetObject();
	  if (object)
		if (!object->Mailed()) {
		  object->Mail();
		  object->mCounter = 1;
		} else
		  object->mCounter++;
	  
#if EVAL_VIEWCELLS	 
	  // and do the same for all viewcells
	  ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
	  
	  for (; it != (*ri).mRay->mViewCells.end(); ++it) {
		ViewCell *viewcell = *it;
		if (!viewcell->Mailed()) {
		  currInfo.viewCells++;
		  viewcell->Mail();
		  viewcell->mCounter = 1;
		} else
		  viewcell->mCounter++;
	  }
#endif
	}

  if (sumWeights!=0.0f)
	currInfo.contribution = sumContribution/sumWeights;
  else
	currInfo.contribution = 0.0f;

  sumWeightsFront = sumWeights;
  sumContributionFront = sumContribution;

  sumWeightsBack = 0;
  sumContributionBack = 0;
  
  currInfo.viewCellsBack = 0;
  currInfo.viewCellsFront = currInfo.viewCells;
  
  Intersectable::NewMail();
  
  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
	  ci < splitCandidates->end();
	  ci++) {
	switch ((*ci).type) {
	case SortableEntry::ERayMin: {
	  currInfo.raysFront--;
	  currInfo.raysBack++;

	  RssTreeNode::RayInfo *rayInfo = (RssTreeNode::RayInfo *) (*ci).data;
	  
	  float weight, contribution;
	  GetRayContribution(*rayInfo, weight, contribution);

	  sumWeightsBack += weight;
	  sumContributionBack += contribution;

	  sumWeightsFront -= weight;
	  sumContributionFront -= contribution;

	  Intersectable *object = rayInfo->GetObject();
	  if (object) {
		if (!object->Mailed()) {
		  object->Mail();
		  currInfo.pvsBack++;
		}
		if (--object->mCounter == 0)
		  currInfo.pvsFront--;
	  }

#if EVAL_VIEWCELLS	 
	  ViewCellContainer::const_iterator it = rayInfo->mRay->mViewCells.begin();
	  for (; it != rayInfo->mRay->mViewCells.end(); ++it) {
		ViewCell *viewcell = *it;
		if (!viewcell->Mailed()) {
		  viewcell->Mail();
		  currInfo.viewCellsBack++;
		}
		if (--viewcell->mCounter == 0)
		  currInfo.viewCellsFront--;
	  }
#endif
	  
	  break;
	}
	}

	float position = (*ci).value;
	
	if (position > minBand && position < maxBand) {
	  currInfo.position = position;

	  
	  if (sumWeightsBack!=0.0f)
		info.contributionBack = sumContributionBack/sumWeightsBack;
	  else
		info.contributionBack = 0.0f;
	  
	  if (sumWeightsFront!=0.0f)
		info.contributionFront = sumContributionFront/sumWeightsFront;
	  else
		info.contributionFront = 0.0f;
	  
	  GetCostRatio(
				   leaf,
				   currInfo);
	  
			
	  //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	  //      cout<<"cost= "<<sum<<endl;
	  
	  if (currInfo.costRatio < info.costRatio) {
		info = currInfo;
      }
    }
  }
  
  
  //  cout<<"===================="<<endl;
  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
}

void
RssTree::SortSubdivisionCandidates(
							 RssTreeLeaf *node,
							 const int axis
							 )
{
  
  splitCandidates->clear();
  
  int requestedSize = 2*(node->rays.size());
  // creates a sorted split candidates array
  if (splitCandidates->capacity() > 500000 &&
      requestedSize < (int)(splitCandidates->capacity()/10) ) {
    
    delete splitCandidates;
    splitCandidates = new vector<SortableEntry>;
  }
  
  splitCandidates->reserve(requestedSize);

  // insert all queries 
  for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
      ri < node->rays.end();
      ri++) {
	if ((*ri).mRay->IsActive()) {
	  if (axis < 3) {
		splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
												 (*ri).GetOrigin(axis),
												 (void *)&(*ri))
								   );
	  } else {
		float pos = (*ri).GetDirParametrization(axis-3);
		splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
												 pos,
												 (void *)&(*ri))
								   );
	  }
	}
  }
  
  stable_sort(splitCandidates->begin(), splitCandidates->end());
}


void
RssTree::EvaluateLeafStats(const TraversalData &data)
{

  // the node became a leaf -> evaluate stats for leafs
  RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;

  if (data.depth >= termMaxDepth)
    stat.maxDepthNodes++;
  
  //  if ( (int)(leaf->rays.size()) < termMinCost)
  //    stat.minCostNodes++;
  if ( leaf->GetPvsSize() <= termMinPvs)
	stat.minPvsNodes++;

  if ( leaf->GetPvsSize() <= termMinRays)
	stat.minRaysNodes++;

  if (leaf->GetAvgRayContribution() > termMaxRayContribution )
	stat.maxRayContribNodes++;
	
  if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
	stat.minSizeNodes++;
  }

  if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
    stat.maxRayRefs = leaf->rays.size();

}

bool
RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
{
  return ( (leaf->GetPvsSize() <= termMinPvs) ||
		   (leaf->rays.size() <= termMinRays) ||
		   //			 (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
		   (leaf->depth >= termMaxDepth) ||
		   (SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize) 
		   );
}


RssTreeNode *
RssTree::SubdivideNode(
					   RssTreeLeaf *leaf,
					   const AxisAlignedBox3 &box,
					   AxisAlignedBox3 &backBBox,
					   AxisAlignedBox3 &frontBBox
					   )
{
  
  if (TerminationCriteriaSatisfied(leaf)) {
#if 0
	if (leaf->depth >= termMaxDepth) {
	  cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
	  cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
	}
#endif
		
	return leaf;
  }
	
  SplitInfo info;
	
  // select subdivision axis
  int axis = SelectPlane( leaf,
						  box,
						  info
						  );
  //  Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
  
  if (axis == -1) {
    return leaf;
  }
  
  stat.nodes+=2;
  stat.leaves += 1;
  stat.splits[axis]++;

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

  node->axis = axis;
  node->position = info.position;

  node->bbox = GetBBox(leaf);
  node->dirBBox = GetDirBBox(leaf);
  
  RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
  RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);

  
  
  // replace a link from node's parent
  if (  leaf->parent )
    leaf->parent->ReplaceChildLink(leaf, node);
  // and setup child links
  node->SetupChildLinks(back, front);

  back->bbox = leaf->bbox;
  front->bbox = leaf->bbox;
  back->dirBBox = leaf->dirBBox;
  front->dirBBox = leaf->dirBBox;
  
  if (axis <= RssTreeNode::SPLIT_Z) {
	back->bbox.SetMax(axis, info.position);
    front->bbox.SetMin(axis, info.position);
  } else {
	back->dirBBox.SetMax(axis-3, info.position);
    front->dirBBox.SetMin(axis-3, info.position);
  }
  
  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	  ri != leaf->rays.end();
	  ri++) {
	if ((*ri).mRay->IsActive()) {
	  
	  // first unref ray from the former leaf
	  (*ri).mRay->Unref();
	  
	  // Debug << "computed t: " << (*ri).mRay->mT << endl;
	  // determine the side of this ray with respect to the plane
	  int side = node->ComputeRaySide(*ri);
	  
	  if (side == 1) 
		front->AddRay(*ri);
	  else
		back->AddRay(*ri);
	} else
	  (*ri).mRay->Unref();
  }

  // distribute the total number of rays according to the distribution
  // of rays which remained
  //  front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
  //  back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
  
#if 0
  front->SetPvsSize(pvsFront);
  back->SetPvsSize(pvsBack);
  // compute entropy as well
  front->ComputeEntropyImportance();
  back->ComputeEntropyImportance();
#else
  UpdatePvsSize(front);
  UpdatePvsSize(back);
#endif

  
  // update stats
  stat.rayRefs -= (int)leaf->rays.size();
  stat.rayRefs += back->rays.size() + front->rays.size();

  
  delete leaf;
  return node;
}






int
RssTree::ReleaseMemory(const int time)
{
  stack<RssTreeNode *> tstack;
  
  // find a node in the tree which subtree will be collapsed
  int maxAccessTime = time - accessTimeThreshold;
  int released;

  PushRoots(tstack);

  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
    
 
    if (!node->IsLeaf()) {
      RssTreeInterior *in = (RssTreeInterior *)node;
      //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
      if (in->depth >= minCollapseDepth &&
		  in->lastAccessTime <= maxAccessTime) {
		released = CollapseSubtree(node, time);
		break;
      }
      
      if (in->back->GetAccessTime() <  
		  in->front->GetAccessTime()) {
		tstack.push(in->front);
		tstack.push(in->back);
      } else {
		tstack.push(in->back);
		tstack.push(in->front);
      }
    }
  }

  while (tstack.empty()) {
    // could find node to collaps...
    //    cout<<"Could not find a node to release "<<endl;
    break;
  }
  
  return released;
}




RssTreeNode *
RssTree::SubdivideLeaf(
					   RssTreeLeaf *leaf
					   )
{
  RssTreeNode *node = leaf;
	
  AxisAlignedBox3 leafBBox = GetBBox(leaf);

  static int pass = 0;
  pass ++;
	
  // check if we should perform a dynamic subdivision of the leaf
  if (!TerminationCriteriaSatisfied(leaf)) {
    
	// memory check and realese...
    if (GetMemUsage() > maxTotalMemory) {
      ReleaseMemory( pass );
    }
    
    AxisAlignedBox3 backBBox, frontBBox;

    // subdivide the node
    node =
      SubdivideNode(leaf,
					leafBBox,
					backBBox,
					frontBBox
					);
  }
	
  return node;
}



void
RssTree::UpdateRays(
					VssRayContainer &remove,
					VssRayContainer &add
					)
{
  RssTreeLeaf::NewMail();

  // schedule rays for removal
  for(VssRayContainer::const_iterator ri = remove.begin();
      ri != remove.end();
      ri++) {
    (*ri)->ScheduleForRemoval();
  }

  int inactive=0;
  Debug<<"Deleting rays..."<<endl<<flush;

  for(VssRayContainer::const_iterator ri = remove.begin();
      ri != remove.end();
      ri++) {
    if ((*ri)->ScheduledForRemoval())
	  //      RemoveRay(*ri, NULL, false);
	  // !!! BUG - with true it does not work correctly - aggreated delete
      RemoveRay(*ri, NULL, true);
    else
      inactive++;
  }

  Debug<<"done."<<endl<<flush;

  //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
  Debug<<"Adding rays..."<<endl<<flush;

  for(VssRayContainer::const_iterator ri = add.begin();
      ri != add.end();
      ri++) {
	RssTreeNode::RayInfo info(*ri);
	//	if (mForcedBoundingBox==NULL || ClipRay(info, bbox))
	AddRay(info);
  }

  Debug<<"done."<<endl<<flush;

  stat.rayRefs += add.size() - remove.size();

  Debug<<"Updating statistics..."<<endl<<flush;
  UpdateTreeStatistics();
  Debug<<"done."<<endl<<flush;

  // check whether the tree should be prunned
  if (stat.rayRefs > mMaxRays) {
	Debug<<"Prunning rays..."<<endl<<flush;
	PruneRays((int)(0.8f*mMaxRays));
	Debug<<"done."<<endl<<flush;
	//	UpdateTreeStatistics();
  }
}

  


void
RssTree::RemoveRay(VssRay *ray,
				   vector<RssTreeLeaf *> *affectedLeaves,
				   const bool removeAllScheduledRays
				   )
{
	
  stack<RayTraversalData> tstack;

  PushRoots(tstack, RssTreeLeaf::RayInfo(ray));
  
  RayTraversalData data;

  // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;

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

    if (!data.node->IsLeaf()) {
      // split the set of rays in two groups intersecting the
      // two subtrees

      TraverseInternalNode(data, tstack);
      
    } else {
      // remove the ray from the leaf
      // find the ray in the leaf and swap it with the last ray...
      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
      
      if (!leaf->Mailed()) {
		leaf->Mail();
		if (affectedLeaves)
		  affectedLeaves->push_back(leaf);
	
		if (removeAllScheduledRays) {
		  int tail = (int)leaf->rays.size()-1;

		  for (int i=0; i < (int)(leaf->rays.size()); i++) {
			if (leaf->rays[i].mRay->ScheduledForRemoval()) {
			  // find a ray to replace it with
			  while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
				stat.removedRayRefs++;
				leaf->rays[tail].mRay->Unref();
				leaf->rays.pop_back();
				tail--;
			  }

			  if (tail < i)
				break;
	      
			  stat.removedRayRefs++;
			  leaf->rays[i].mRay->Unref();
			  leaf->rays[i] = leaf->rays[tail];
			  leaf->rays.pop_back();
			  tail--;
			}
		  }
		}
      }
      
      if (!removeAllScheduledRays)
		for (int i=0; i < (int)leaf->rays.size(); i++) {
		  if (leaf->rays[i].mRay == ray) {
			stat.removedRayRefs++;
			ray->Unref();
			leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
			leaf->rays.pop_back();
			// check this ray again
			break;
		  }
		}
      
    }
  }

  if (ray->RefCount() != 0) {
    cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
    exit(1);
  }
  
}


void
RssTree::AddRay(RssTreeNode::RayInfo &info)
{
  if (EpsilonEqualV3(info.GetOrigin(), info.GetTermination(), Limits::Small))
	return;

  stack<RayTraversalData> tstack;

  RssTreeNode *root = GetRoot(info.GetSourceObject());
  tstack.push(RayTraversalData(root, info));
  
  RayTraversalData data;
  

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

    if (!data.node->IsLeaf()) {
      TraverseInternalNode(data, tstack);
    } else {
      // remove the ray from the leaf
      // find the ray in the leaf and swap it with the last ray...
      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
      leaf->AddRay(data.rayData);
      stat.addedRayRefs++;
    }
  }
}

void
RssTree::TraverseInternalNode(
							  RayTraversalData &data,
							  stack<RayTraversalData> &tstack)
{
  RssTreeInterior *in =  (RssTreeInterior *) data.node;
  
  // determine the side of this ray with respect to the plane
  int side = in->ComputeRaySide(data.rayData
								);
  
  if (side == 1)
	tstack.push(RayTraversalData(in->front, data.rayData));
  else
	tstack.push(RayTraversalData(in->back, data.rayData));
  
}


int
RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
{
  // first count all rays in the subtree
  // use mail 1 for this purpose
  stack<RssTreeNode *> tstack;
  int rayCount = 0;
  int totalRayCount = 0;
  int collapsedNodes = 0;

#if DEBUG_COLLAPSE
  cout<<"Collapsing subtree"<<endl;
  cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
  cout<<"depth="<<(int)sroot->depth<<endl;
#endif

  //    tstat.collapsedSubtrees++;
  //    tstat.collapseDepths += (int)sroot->depth;
  //    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
  
  tstack.push(sroot);
  VssRay::NewMail();
	
  while (!tstack.empty()) {
    collapsedNodes++;
    RssTreeNode *node = tstack.top();
    tstack.pop();

    if (node->IsLeaf()) {
      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		  ri != leaf->rays.end();
		  ri++) {
				
		totalRayCount++;
		if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
		  (*ri).mRay->Mail();
		  rayCount++;
		}
      }
    } else {
      tstack.push(((RssTreeInterior *)node)->back);
      tstack.push(((RssTreeInterior *)node)->front);
    }
  }
  
  VssRay::NewMail();

  // create a new node that will hold the rays
  RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
  if (  newLeaf->parent )
    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
  

  tstack.push( sroot );

  while (!tstack.empty()) {

    RssTreeNode *node = tstack.top();
    tstack.pop();

    if (node->IsLeaf()) {
      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
      
      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		  ri != leaf->rays.end();
		  ri++) {
				
		// unref this ray from the old node
				
		if ((*ri).mRay->IsActive()) {
		  (*ri).mRay->Unref();
		  if (!(*ri).mRay->Mailed()) {
			(*ri).mRay->Mail();
			newLeaf->AddRay(*ri);
		  }
		} else
		  (*ri).mRay->Unref();
				
      }
    } else {
      tstack.push(((RssTreeInterior *)node)->back);
      tstack.push(((RssTreeInterior *)node)->front);
    }
  }
  
  // delete the node and all its children
  delete sroot;
  
  //   for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
  //       ri != newLeaf->rays.end();
  //       ri++)
  //     (*ri).ray->UnMail(2);


#if DEBUG_COLLAPSE
  cout<<"Total memory before="<<GetMemUsage()<<endl;
#endif

  stat.nodes -= collapsedNodes - 1;
  stat.leaves -= collapsedNodes/2 - 1;
  stat.rayRefs -= totalRayCount - rayCount;
  
#if DEBUG_COLLAPSE
  cout<<"collapsed nodes"<<collapsedNodes<<endl;
  cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
  cout<<"Total memory after="<<GetMemUsage()<<endl;
  cout<<"================================"<<endl;
#endif

  //  tstat.collapsedNodes += collapsedNodes;
  //  tstat.collapsedRays += totalRayCount - rayCount;
    
  return totalRayCount - rayCount;
}


int
RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
{
  stack<RssTreeNode *> tstack;
  PushRoots(tstack);

  Intersectable::NewMail();
  int pvsSize = 0;
	
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
    
 
    if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		  ri != leaf->rays.end();
		  ri++)
		if ((*ri).mRay->IsActive()) {
		  Intersectable *object;
		  object = (*ri).GetObject();
		  if (object && !object->Mailed()) {
			pvsSize++;
			object->Mail();
		  }
		}
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  if (in->axis < 3) {
		if (box.Max(in->axis) >= in->position )
		  tstack.push(in->front);
				
		if (box.Min(in->axis) <= in->position )
		  tstack.push(in->back);
	  } else {
		// both nodes for directional splits
		tstack.push(in->front);
		tstack.push(in->back);
	  }
	}
  }
  return pvsSize;
}

int
RssTree::CollectPvs(const AxisAlignedBox3 &box,
					ObjectContainer &pvs
					) const
{
  stack<RssTreeNode *> tstack;
  PushRoots(tstack);
  
  Intersectable::NewMail();
  
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
    
 
    if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		  ri != leaf->rays.end();
		  ri++)
		if ((*ri).mRay->IsActive()) {
		  Intersectable *object;
		  object = (*ri).GetObject();
		  if (object && !object->Mailed()) {
			pvs.push_back(object);
			object->Mail();
		  }
		}
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  if (in->axis < 3) {
		if (box.Max(in->axis) >= in->position )
		  tstack.push(in->front);
		
		if (box.Min(in->axis) <= in->position )
		  tstack.push(in->back);
	  } else {
		// both nodes for directional splits
		tstack.push(in->front);
		tstack.push(in->back);
	  }
	}
  }
  return pvs.size();
}

void
RssTree::UpdateTreeStatistics()
{
  for (int i=0; i < mRoots.size(); i++)
	UpdateImportance(mRoots[i]);
  
#if 0  
  stack<RssTreeNode *> tstack;
	
  float sumPvsSize = 0.0f;
  float sumRayContribution = 0.0f;
  float sumWeightedRayContribution = 0.0f;
  float sumImportance = 0.0f;
  float sumPvsEntropy = 0.0f;
  float sumRayLengthEntropy = 0.0f;
  float sumRays = 0.0f;
  
  int leaves = 0;


  PushRoots(tstack);
  
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
	tstack.pop();
	
    if (node->IsLeaf()) {
	  leaves++;
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  // updated above already
	  //	  UpdatePvsSize(leaf);
	  
	  sumPvsSize += leaf->GetPvsSize();
	  sumRayContribution += leaf->GetAvgRayContribution();
	  

	  RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
	  for (;it != leaf->rays.end(); ++it) {
		float weight, contribution;
		GetRayContribution(*it, weight, contribution);
		sumWeightedRayContribution += weight*contribution;
	  }
	  
	  //	  sumPvsEntropy += leaf->mPvsEntropy;
	  //	  sumRayLengthEntropy += leaf->mRayLengthEntropy;
	  sumRays += leaf->rays.size();
	  
	  float imp = leaf->GetImportance();
	  
	  //	  if (imp > 1.0f)
	  //		cout<<"warning imp > 1.0f:"<<imp<<endl;
	  
	  sumImportance += imp;
	  
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }
  
  stat.avgPvsSize = sumPvsSize/(float)leaves;
  stat.avgRays = sumRays/(float)leaves;
  stat.avgRayContribution = sumRayContribution/(float)leaves;
  //  avgPvsEntropy = sumPvsEntropy/(float)leaves;
  //  avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
  stat.avgImportance = sumImportance/(float)leaves;
  stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
  stat.rayRefs = (int)sumRays;
#endif
  
}

void
RssTree::UpdateImportance(
						  RssTreeNode *node)
{
  if (node->IsLeaf()) {
	UpdatePvsSize((RssTreeLeaf *)node);
  } else {
	RssTreeInterior *in = (RssTreeInterior *)node;
	UpdateImportance(in->front);
	UpdateImportance(in->back);
	node->mImportance = in->front->mImportance + in->back->mImportance;
  }
}

// generate a single ray described by parameters in a unit cube
void
RssTree::GenerateRay(float *params, SimpleRayContainer &rays)
{
  RssTreeNode *node;
  
  // works only with single root now!
  node = mRoots[0];
  while (!node->IsLeaf()) {
	RssTreeInterior *in = (RssTreeInterior *)node;
	// decide whether to go left or right
	int axis = in->axis;
	//	float p = in->GetRelativePosition(); 
	float p = in->back->GetImportance()/(in->back->GetImportance() + in->front->GetImportance());

	//	cout<<"p = "<<p<<endl;
	//	cout<<"axis="<<axis<<" pos="<<pos<<endl;
	//	cout<<GetBBox(in)<<endl;
	//	cout<<GetDirBBox(in)<<endl;
	//	cout<<"********"<<endl;
	
	if (params[axis] < p) {
	  node = in->back;
	  if (p > Limits::Small)
		params[axis] = params[axis] / p;
	} else {
	  node = in->front;
	  if (p < 1.0f - Limits::Small)
		params[axis] = (params[axis] - p) / (1.0f - p);
	}
  }

  //  cout<<params[0]<<" "<<params[1]<<" "<<params[2]<<" "<<params[3]<<" "<<params[4]<<endl;
  
  GenerateLeafRay(params, rays, node->bbox, node->dirBBox);
}

int
RssTree::GenerateRays(
					  const float ratioPerLeaf,
					  SimpleRayContainer &rays)
{

  stack<RssTreeNode *> tstack;
  PushRoots(tstack);
	
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
		
    if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  float c = leaf->GetImportance();
	  int num = int(c*ratioPerLeaf + 0.5f);
	  //			cout<<num<<" ";

	  for (int i=0; i < num; i++) {
		Vector3 origin = GetBBox(leaf).GetRandomPoint();
		Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
		Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
		//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
		rays.push_back(SimpleRay(origin,
								 direction,
								 SamplingStrategy::RSS_BASED_DISTRIBUTION,
								 1.0f
								 ));
	  }
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }

  return (int)rays.size();
}

void
RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
{
  stack<RssTreeNode *> tstack;
  PushRoots(tstack);
	
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
		
    if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  leaves.push_back(leaf);
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }
}

bool
RssTree::ValidLeaf(RssTreeLeaf *leaf) const
{
  return true;
  //return leaf->rays.size() > termMinRays/4;
}


void
RssTree::PickEdgeRays(RssTreeLeaf *leaf,
					  int &r1,
					  int &r2
					  )
{
  int nrays = leaf->rays.size();

  if (nrays == 2) {
	r1 = 0;
	r2 = 1;
	return;
  }

#if 1
  int tries = min(20, nrays);

  while (--tries) {
	r1 = Random(nrays);
	r2 = Random(nrays);
	if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
	  break;
  }

  if (r1 == r2)
	r2 = (r1+1)%leaf->rays.size();

#else
  // pick a random ray
  int base = Random(nrays);
  RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
  Intersectable *baseObject = baseRay->GetObject();
  
  // and a random 5D derivative which will be used to find the minimal projected distances
  Vector3 spatialDerivative;
  Vector3 directionalDerivative;

  while (1) {
	spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
								RandomValue(-1.0f, 1.0f),
								RandomValue(-1.0f, 1.0f));
	float m = Magnitude(spatialDerivative);
	if (m != 0) {
	  spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
	  break;
	}
  }

  while (1) {
	directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
									RandomValue(-1.0f, 1.0f),
									0.0f);
	float m = Magnitude(directionalDerivative);
	if (m != 0) {
	  directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
	  break;
	}
  }

  // now find the furthest sample from the same object and the closest from a different object
  int i = 0;
  float minDist = MAX_FLOAT;
  float maxDist = -MAX_FLOAT;
  r1 = base;
  r2 = base;
  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	  ri != leaf->rays.end();
	  ri++, i++) {
	float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
													 *ri,
													 spatialDerivative,
													 directionalDerivative
													 );

	if ((*ri).GetObject() == baseObject) {
	  if (dist > maxDist) {
		maxDist = dist;
		r1 = i;
	  }
	} else {
	  if (dist > 0.0f && dist < minDist) {
		minDist = dist;
		r2 = i;
	  }
	}
  }
  
#endif
}


struct RayNeighbor {
  int rayInfo;
  float distance;
  RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
};

void
RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
							vector<RssTreeLeaf::SilhouetteRays> &rays
							)
{
  // for every leaf find its neares neighbor from a different object
  vector<RayNeighbor> neighbors(leaf->rays.size());
  int i, j;
  for (i=0; i < leaf->rays.size(); i++)
	for (j=0; j < leaf->rays.size(); j++)
	  if (i!=j) {
		float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
		if (d < neighbors[i].distance) {
		  neighbors[i].distance = d;
		  neighbors[i].rayInfo = j;
		}
	  }

  // now check which are the pairs of nearest neighbors
  for (i=0; i < leaf->rays.size(); i++) {
	int j = neighbors[i].rayInfo;
	if (neighbors[j].rayInfo == i) {
	  // this is a silhouette edge pair
	  if (i < j) {
		// generate an silhoutte ray pair
		rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
												   &leaf->rays[j]));
	  }
	} else {
	  // this is not a silhouette ray - delete???
	}
  }
  
}


bool
GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
				 const int index,
				 int &i1,
				 int &i2,
				 int &i3)
{
  int found = 0;
  int indices[3];

  int size = rays.size();
  // use russian roulete selection for the tripple
  // number of free positions for the bullet
  int positions = size - 1;
  int i;
  for (i=0; i < size; i++) {
	if (rays[i]->mRay->Mailed())
	  positions--;
  }
  
  if (positions < 3)
	return false;
  
  for (i=0; i < size; i++) {
	if (i != index && !rays[i]->mRay->Mailed()) {
	  float p = (3 - found)/(float)positions;
	  if (Random(1.0f) < p) {
		indices[found] = i;
		found++;
	  }
	}
	positions--;
  }
  return true; // corr. matt
}

bool
RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
								const RssTreeNode::RayInfo &r1,
								const RssTreeNode::RayInfo &r2,
								const RssTreeNode::RayInfo &r3)
{
  

  return false;
}

int
RssTree::RemoveInteriorRays(
							RssTreeLeaf *leaf
							)
{
#if 1
  // first collect all objects refered in this leaf
  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> > rayMap;
  
  RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
  for (; it != leaf->rays.end(); ++it) {
	Intersectable *object = (*it).GetObject();
	
	rayMap[object].push_back(&(*it));
// 	vector<RayInfo *> *data = rayMap.Find(object);
// 	if (data) {
// 	  data->push_back(&(*it));
// 	} else {
// 	  //	  rayMap[object] = vector<RayInfo *>;
// 	  rayMap[object].push_back(&(*it));
// 	}
  }

  // now go through all objects
  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> >::iterator mi;

  // infos of mailed rays are scheduled for removal
  VssRay::NewMail();
  for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
	vector<RssTreeLeaf::RayInfo *> &rays = (*mi).second;
	vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
	int rayIndex = 0;
	
	for (; ri != rays.end(); ++ri, ++rayIndex) {
	  RssTreeNode::RayInfo *ray = *ri;
	  int tries = rays.size();
	  for (int i = 0; i < tries; i++) {
		int r1, r2, r3;
		GetRandomTripple(rays,
						 rayIndex,
						 r1,
						 r2,
						 r3);
		if (IsRayConvexCombination(*ray,
								   *rays[r1],
								   *rays[r2],
								   *rays[r3])) {
		  ray->mRay->Mail();
		}
	  }
	}
  }
  

#endif  
	return 0;
}


void
RssTree::GenerateLeafRaysSimple(RssTreeLeaf *leaf,
								const int numberOfRays,
								SimpleRayContainer &rays)
{
  
  if (numberOfRays == 0)
	return;
  
  int startIndex = (int)rays.size();
  //  Debug<<"B"<<flush;
  
  AxisAlignedBox3 box = GetBBox(leaf);
  AxisAlignedBox3 dirBox = GetDirBBox(leaf);

  float boxSize = Magnitude(box.Diagonal());
  
  const float smoothRange = 0.1f;
  if (smoothRange != 0.0f) {
	box.Scale(1.0f + smoothRange);
	dirBox.Scale(1.0f + smoothRange);
  }

  //  Debug<<"R"<<flush;
  
  for (int i=0; i < numberOfRays; i++) {
	//	Debug<<i<<flush;
	float r[5];
	
	r[0] = RandomValue(0.0f, 1.0f);
	r[1] = RandomValue(0.0f, 1.0f);
	r[2] = RandomValue(0.0f, 1.0f);
	r[3] = RandomValue(0.0f, 1.0f);
	r[4] = RandomValue(0.0f, 1.0f);
	
	GenerateLeafRay(r, rays, box, dirBox);
  }

  //  Debug<<"D"<<flush;
  
  float density = 1.0f;
  float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
  for (int i=startIndex; i < rays.size(); i++) {
	rays[i].mPdf = p;
  }
}

void
RssTree::GenerateLeafRay(
						 float *params,
						 SimpleRayContainer &rays,
						 const AxisAlignedBox3 &sbox,
						 const AxisAlignedBox3 &sDirBBox
						 )
{
  Vector3 origin = sbox.GetPoint(Vector3(params[0], params[1], params[2]));
  Vector3 dirVector = sDirBBox.GetPoint(Vector3(params[3], params[4], 0.0f));
  Vector3 direction = VssRay::GetInvDirParam(dirVector.x, dirVector.y);
  
  //	float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
  float dist = 0.0f;
  float minT, maxT;
  
  // compute interection of the ray with the box
#if 1
  if (sbox.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
	dist = (maxT + Limits::Small);
#else
  float boxSize = Magnitude(sbox.Size());
  dist = 0.5f*boxSize;
#endif
  
  origin += direction*dist;
  
  //Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
  rays.push_back(SimpleRay(origin,
						   direction,
						   SamplingStrategy::RSS_BASED_DISTRIBUTION,
						   1.0f
						   ));
}

void
RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
						  const int numberOfRays,
						  SimpleRayContainer &rays)
{

  if (numberOfRays == 0)
	return;
  
  
  int nrays = (int)leaf->rays.size();
  int startIndex = (int)rays.size();
  
  AxisAlignedBox3 box = GetBBox(leaf);
  AxisAlignedBox3 dirBox = GetDirBBox(leaf);


  const float smoothRange = 0.1f;
  if (smoothRange != 0.0f) {
	box.Scale(1.0f + smoothRange);
	dirBox.Scale(1.0f + smoothRange);
  }
  
  Exporter *exporter = NULL;
  VssRayContainer selectedRays;
  static int counter = 0;
  if (counter < 0) {
	char filename[128];
	sprintf(filename, "selected-rays%04d.x3d", counter);
	exporter = Exporter::GetExporter(filename);
	//	exporter->SetWireframe();
	//	exporter->ExportKdTree(*mKdTree);
	exporter->SetWireframe();
//	exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
	exporter->SetWireframe();
	exporter->ExportBox(box);
	exporter->SetFilled();
	counter++;
  }
  

  float extendedConvexCombinationProb = 0.0f; //0.7f
  //  float silhouetteCheckPercentage = 1.0f; //0.5f
  float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
  float silhouetteObjectCheckPercentage = 0.8f;

  //  int numberOfTries = numberOfRays*4;
  int numberOfTries = numberOfRays;
  int generated = 0;

  vector<Intersectable *> pvs(0);
  if (silhouetteCheckPercentage != 0.0f) {
	Intersectable::NewMail();
	pvs.reserve(leaf->GetPvsSize());
	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		ri != leaf->rays.end();
		ri++) {
	  Intersectable *object = (*ri).GetObject();
	  if (object) {
		if (!object->Mailed()) {
		  object->Mail();
		  object->mCounter = 1;
		  pvs.push_back(object);
		  if (exporter)
			exporter->ExportIntersectable(object);
		} else {
		  object->mCounter++;
		}
	  }
	}
	// sort objects based on their mCounter
	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
  }
  
  for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
														  extendedConvexCombinationProb));
	
	Vector3 origin, direction;
	// generate half of convex combination and half of random rays
	if (useExtendedConvexCombination) {
	  // pickup 2 random rays
	  int r1, r2;
	  PickEdgeRays(leaf, r1, r2);

	  
	  Vector3 o1 = leaf->rays[r1].GetOrigin();
	  Vector3 o2 = leaf->rays[r2].GetOrigin();
	  
	  const float overlap = 0.0f;
	  float w1, w2;
	  GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
	  origin = w1*o1 + w2*o2;
	  GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
	  direction =
		w1*leaf->rays[r1].GetDir() +
		w2*leaf->rays[r2].GetDir();

	  // shift the origin a little bit
	  direction.Normalize();
	  
	} else {
	  Vector3 dirVector;

	  Vector3 pVector;
	  Vector3 dVector;

	  pVector = Vector3(RandomValue(0.0f, 1.0f),
						RandomValue(0.0f, 1.0f),
						RandomValue(0.0f, 1.0f));
	  
	  dVector = Vector3(RandomValue(0.0f, 1.0f),
						RandomValue(0.0f, 1.0f),
						0.0f);
	  
	  
	  origin = box.GetPoint(pVector);
	  dirVector = dirBox.GetPoint(dVector);

	  direction = VssRay::GetInvDirParam(dirVector.x,dirVector.y);
	}
	
	
	//	float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
	float dist = 0.0f;
	float minT, maxT;
	// compute interection of the ray with the box

	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
	  dist = maxT;
  
	
	origin += direction*dist;
	
	bool intersects = false;
	
	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
	  if (exporter) {
		VssRay *ray = new VssRay(origin, origin + 100*direction, NULL, NULL);
		selectedRays.push_back(ray);
	  }
	  
	  // check whether the ray does not intersect already visible objects
	  Ray traversalRay;
	  traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
	  for(vector<Intersectable *>::const_iterator oi = pvs.begin();
		  oi != pvs.end();
		  oi++) {
		Intersectable *object = *oi;
		// do not test every object
		if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
			 object->CastRay(traversalRay) ) {
		  intersects = true;
		  break;
		}
	  }
	}
	
	
	if (!intersects) {
	  //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	  rays.push_back(SimpleRay(origin,
							   direction,
							   SamplingStrategy::RSS_BASED_DISTRIBUTION,
							   1.0f
							   ));
	  generated++;
	}
  }

  //  cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
  // assign pdfs to the generated rays
  // float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
  float density = 1.0f;

  float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
  for (int i=startIndex; i < rays.size(); i++) {
	rays[i].mPdf = p;
  }
  
  
  if (exporter) {
	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
	delete exporter;
	CLEAR_CONTAINER(selectedRays);
  }
  
}

int
RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
						 const float ratio)
{
  int i;
  int j;
  
  for (j=0, i=0; i < leaf->rays.size(); i++) {
	if (Random(1.0f) < ratio) {
	  // copy a valid sample
	  if (i!=j) 
		leaf->rays[j] = leaf->rays[i];
	  j++;
	} else {
	  // delete the ray
	  leaf->rays[i].mRay->Unref();
	  if (leaf->rays[i].mRay->RefCount() != 0) {
		cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
		exit(1);
	  }
	  delete leaf->rays[i].mRay;
	}
  }


  leaf->rays.resize(j);
  int removed = (i-j);
  stat.rayRefs -= removed;
  return removed;
}

int
RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
							   const float ratio)
{
  int i;
  
  if (leaf->rays.size() == 0)
	return 0;
  
  sort(leaf->rays.begin(),
	   leaf->rays.end(),
	   RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);

  int desired = int(ratio*leaf->rays.size());
  int removed = (int)leaf->rays.size() - desired;
  
  for (i=desired; i < (int)leaf->rays.size(); i++) {
	// delete the ray
	leaf->rays[i].mRay->Unref();
	if (leaf->rays[i].mRay->RefCount() != 0) {
	  cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
	  cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
	  exit(1);
	}
	delete leaf->rays[i].mRay;
  }
  
  leaf->rays.resize(desired);
  stat.rayRefs -= removed;
  return removed;
}


int
RssTree::PruneRays(
				   const int maxRays
				   )
{

  stack<RssTreeNode *> tstack;
  int prunned = 0;

  // prune more rays to amortize  the procedure
  int desired = (int)(maxRays * 0.8f);

  Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;

  // prune random rays from each leaf so that the ratio's remain the same
  PushRoots(tstack);
  
  if ( maxRays >= stat.rayRefs )
	return 0;
  
  float ratio = desired/(float)stat.rayRefs;
  
  while (!tstack.empty()) {
	RssTreeNode *node = tstack.top();
	tstack.pop();
	
	if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  //prunned += PruneRaysRandom(leaf, ratio);
	  prunned += PruneRaysContribution(leaf, ratio);
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }
  
  Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
  
  return prunned;
}

int
RssTree::GenerateRays(const int numberOfRays,
					  const int numberOfLeaves,
					  SimpleRayContainer &rays)
{

  float param[5];
  
  for (int i=0; i < numberOfRays; i++) {
	mHalton.GetNext(5, param);
	GenerateRay(param, rays);
  }
  
  return rays.size();
  
}


float
RssTree::GetAvgPvsSize()
{
  stack<RssTreeNode *> tstack;
  PushRoots(tstack);

  int sumPvs = 0;
  int leaves = 0;
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
		
    if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  // update pvs size
	  UpdatePvsSize(leaf);
	  sumPvs += leaf->GetPvsSize();
	  leaves++;
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }


  return sumPvs/(float)leaves;
}

float weightAbsContributions = ABS_CONTRIBUTION_WEIGHT;
// if small very high importance of the last sample
// if 1.0f then weighs = 1 1/2 1/3 1/4
float passSampleWeightDecay = 2.0f;
//float passSampleWeightDecay = 1000.0f;

float
RssTree::GetSampleWeight(const int pass)
{
  int passDiff = mCurrentPass - pass;
  if (passDiff < 0) {
	cerr<<"Error: passdif < 0"<<passDiff<<endl;
	exit(1);
  }
  float weight;
  if (1) 
	weight = 1.0f/(passDiff + passSampleWeightDecay);
  else 
	switch (passDiff) {
	case 0:
	  //	  weight = 1.0f;
	  //	  break;
	default:
	  weight = 1.0f;
	  break;
	  // 	case 1:
	  // 	  weight = 0.5f;
	  // 	  break;
	  // 	case 2:
	  // 	  weight = 0.25f;
	  // 	  break;
	  // 	case 3:
	  // 	  weight = 0.12f;
	  // 	  break;
	  // 	case 4:
	  // 	  weight = 0.06f;
	  // 	  break;
	  // 	default:
	  // 	  weight = 0.03f;
	  // 	  break;
	}
  return weight;
}

void
RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
							float &weight,
							float &contribution)
{
  VssRay *ray = info.mRay;
  weight = GetSampleWeight(mCurrentPass);
  contribution =
	weightAbsContributions*ray->mPvsContribution +
	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
  // store the computed value
  info.mRay->mWeightedPvsContribution = weight*contribution;
}


void
RssTree::ComputeImportance(RssTreeLeaf *leaf) 
{
#if 0
  if (leaf->mTotalRays)
	leaf->mImportance = leaf->mTotalContribution / leaf->mTotalRays;
  else
	leaf->mImportance = Limits::Small;
  return;
#endif
  
  if (0) 
	leaf->mImportance = leaf->GetAvgRayContribution();
  else {
	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
	  it_end = leaf->rays.end();
	
	float sumContributions = 0.0f;
	float sumRelContributions = 0.0f;
	float sumWeights = 0.0f;
	float maxContribution = 0.0f;
	float maxRelContribution = 0.0f;
	float sumLength = 0;
	
	for (; it != it_end; ++it) {
	  VssRay *ray = (*it).mRay;
	  float sumLength = ray->Length();

	  float weight = GetSampleWeight(ray->mPass);
	  
	  float c1 = weight*ray->mPvsContribution;
	  float c2 = weight*ray->mRelativePvsContribution;

	  sumContributions += c1;
	  if (c1 > maxContribution)
		maxContribution = c1;
	  
	  //$$ 20.7. changed to sqr to pronouce higher contribution so that
	  // they do not get smoothed!!
	  //sumRelContributions += weight*ray->mRelativePvsContribution;
	  
	  sumRelContributions += c2;
	  if (c2 > maxRelContribution)
		maxRelContribution = c2;
	  
	  //sumWeights += weight;
	  sumWeights += 1.0f;
	}

	float distImportance = 0.0f;

	float spatialArea = GetBBox(leaf).SurfaceArea();
	float dirArea = GetDirBBox(leaf).SurfaceArea();
	if (leaf->rays.size()) {
	  // compute average length of the ray 
	  float avgLength = sumLength/leaf->rays.size();
	  // compute estimated area of the visible surface corresponding to these rays
	  float ss = spatialArea/2.0f;
	  float sd = dirArea;
	  float s = ss + avgLength*(2*sqrt(ss*sd) + avgLength*sd);
	  
	  distImportance = s/leaf->rays.size();
	}
	
	// $$ 
	// sumWeights = leaf->mTotalRays;

	//$$ test 30.11. 2006 
	// normalize per volume of the cell not the number of rays!
	// sumWeights = spatialArea*dirArea;
	
	if (sumWeights != 0.0f) {
#if 1
	  leaf->mImportance =
		pow(((weightAbsContributions*sumContributions +
			  (1.0f - weightAbsContributions)*sumRelContributions)/sumWeights), 1.0f);
#else 
	  leaf->mImportance =
		(weightAbsContributions*maxContribution +
		 (1.0f - weightAbsContributions)*maxRelContribution)/sumWeights;
#endif
	} else
	  leaf->mImportance = 0.0f;
	
  }
  
  // return GetAvgRayContribution()*mImportance;
  //return GetAvgRayContribution();
}


float
RssTreeNode::GetImportance()  const
{
  // return mImportance*mImportance;
  return mImportance;
}


float
RssTreeLeaf::ComputePvsEntropy() const
{
  int samples = 0;
  Intersectable::NewMail();
  // set all object as belonging to the fron pvs
  for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
	  ri != rays.end();
	  ri++)
	if ((*ri).mRay->IsActive()) {
	  Intersectable *object = (*ri).GetObject();
	  if (object) {
		if (!object->Mailed()) {
		  object->Mail();
		  object->mCounter = 1;
		} else
		  object->mCounter++;
		samples++;
	  }
	}
  
  float entropy = 0.0f;
  
  if (samples > 1) {
	Intersectable::NewMail();
	for(RayInfoContainer::const_iterator ri = rays.begin();
		ri != rays.end();
		ri++) 
	  if ((*ri).mRay->IsActive()) {
		Intersectable *object = (*ri).GetObject();
		if (object) {
		  if (!object->Mailed()) {
			object->Mail();
			float p = object->mCounter/(float)samples;
			entropy -= p*log(p);
		  }
		}
	  }
	entropy = entropy/log((float)samples);
  }
  else
	entropy = 1.0f;
  
  return entropy;
}

float
RssTreeLeaf::ComputeRayLengthEntropy() const
{
  // get sum of all ray lengths
  // consider only passing rays or originating rays
  float sum = 0.0f;
  int samples = 0;
  int i=0;
  for(RayInfoContainer::const_iterator ri = rays.begin();
      ri != rays.end();
      ri++) 
	if ((*ri).mRay->IsActive()) {
//  	  float s;
//  	  if (i == 0)
//  		s = 200;
//  	  else
//  		s = 100;
//  	  i++;
	  
	  sum += (*ri).mRay->GetSize();
	  samples++;
	}
  
  
  float entropy = 0.0f;
  
  if (samples > 1) {
	i = 0;
	for(RayInfoContainer::const_iterator ri = rays.begin();
		ri != rays.end();
		ri++)
	  if ((*ri).mRay->IsActive()) {
//  		float s;
//  		if (i==0)
//  		  s = 200;
//  		else
//  		  s = 100;
//  		i++;
//  		float p = s/sum;
		float p = (*ri).mRay->GetSize()/sum;
		entropy -= p*log(p);
	  }
	entropy = entropy/log((float)samples);
  } else
	entropy = 1.0f;
  
  return entropy;
}
  


float
RssTreeLeaf::ComputeEntropyImportance() const
{
  
  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
  return 1.0f - ComputePvsEntropy();
}

float
RssTreeLeaf::ComputeRayContributionImportance() const
{
  //  mPvsEntropy = ComputePvsEntropy();
  //  mRayLengthEntropy = ComputeRayLengthEntropy();
  
  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
  return 1.0f - ComputePvsEntropy();
}


AxisAlignedBox3
RssTree::GetShrankedBBox(const RssTreeNode *node) const
{
  if (node->parent == NULL)
	return bbox;
  
  if (!node->IsLeaf())
	return ((RssTreeInterior *)node)->bbox;

  // evaluate bounding box from the ray origins
  AxisAlignedBox3 box;
  box.Initialize();
  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	  ri != leaf->rays.end();
	  ri++)
	if ((*ri).mRay->IsActive()) {
	  box.Include((*ri).GetOrigin());
	}
  return box;
}


int
RssTree::CollectRays(VssRayContainer &rays,
					 const int number)
{
  VssRayContainer allRays;
  CollectRays(allRays);

  
  int desired = min(number, (int)allRays.size());
  float prob = desired/(float)allRays.size();
  //int skip = allRays.size()/desired;
  while (rays.size() < desired) {
	VssRayContainer::const_iterator it = allRays.begin();
	for (; it != allRays.end() && rays.size() < desired; it++) {
	  if (Random(1.0f) < prob)
		rays.push_back(*it);
	}
  }
  return (int)rays.size();
}


int
RssTree::CollectRays(VssRayContainer &rays
					 )
{
  VssRay::NewMail();

  stack<RssTreeNode *> tstack;
  PushRoots(tstack);
  
  while (!tstack.empty()) {
    RssTreeNode *node = tstack.top();
    tstack.pop();
	
    if (node->IsLeaf()) {
	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
	  // update pvs size
	  RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
	  for (;it != leaf->rays.end(); ++it)
		if (!(*it).mRay->Mailed()) {
		  (*it).mRay->Mail();
		  rays.push_back((*it).mRay);
		}
	} else {
	  RssTreeInterior *in = (RssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }
  
  return (int)rays.size();
}


RssTreeNode *
RssTree::GetRoot(Intersectable *object) const
{
  if (mPerObjectTree && object) {
	int id = object->GetId()-1;
	if (id < 0 || id >= mRoots.size()) {
	  Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
		endl<<flush;
	  id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
	}
	return mRoots[id];
  } else
	return mRoots[0];
}

void
RssTree::PushRoots(priority_queue<TraversalData> &st) const
{
  for (int i=0; i < mRoots.size(); i++)
	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
}

void
RssTree::PushRoots(stack<RssTreeNode *> &st)  const
{
  for (int i=0; i < mRoots.size(); i++)
	st.push(mRoots[i]);
}

void
RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info)  const
{
  for (int i=0; i < mRoots.size(); i++)
	st.push(RayTraversalData(mRoots[i], info));
}



bool
RssBasedDistribution::GenerateSample(SimpleRay &ray)
{
  float r[5];

  mHalton.GetNext(5, r);

  if (mRssTree == NULL) {
	// use direction based distribution
	Vector3 origin, direction;
	mPreprocessor.mViewCellsManager->GetViewPoint(origin, Vector3(r[0], r[1], r[2]));
	
	direction = UniformRandomVector(r[3], r[4]);

	const float pdf = 1.0f;
	ray = SimpleRay(origin, direction, RSS_BASED_DISTRIBUTION, pdf);
	
  } else {
	
	SimpleRayContainer rays;
	mRssTree->GenerateRay(r, rays);
	ray = rays[0];
  }
  
  return true;
}

void
RssBasedDistribution::Update(VssRayContainer &vssRays)
{
  if (mRssTree == NULL) {
	mRssTree = new RssTree;
	mRssTree->SetPass(mPass);
	
	/// compute view cell contribution of rays if view cells manager already constructed
	//  mViewCellsManager->ComputeSampleContributions(mVssRays, true, false);
	mRssTree->Construct(mPreprocessor.mObjects, vssRays);
	mRssTree->stat.Print(mPreprocessor.mStats);
	cout<<"RssTree root PVS size = "<<mRssTree->GetRootPvsSize()<<endl;
	
  } else {
	mRssTree->SetPass(mPass);
	Debug<<"Adding rays...\n"<<flush;
	mRssTree->AddRays(vssRays);
	Debug<<"done.\n"<<flush;
	if (1) {
	  //	if (mUpdateSubdivision) {
	  Debug<<"Updating rss tree...\n"<<flush;
	  int subdivided = mRssTree->UpdateSubdivision();
	  Debug<<"done.\n"<<flush;
	  cout<<"subdivided leafs = "<<subdivided<<endl;
	  cout<<"#total leaves = "<<mRssTree->stat.Leaves()<<endl;
	}
  }
  mPass++;
}

}