// ================================================================
// $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 "VssTree.h"

#include "Environment.h"
#include "VssRay.h"
#include "Intersectable.h"
#include "Ray.h"

namespace GtpVisibilityPreprocessor {

#define DEBUG_SPLIT_COST 0

// Static variables
int
VssTreeLeaf::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);
	}
  }
}

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

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


  float refDirAngle;
  Environment::GetSingleton()->GetFloatValue("VssTree.refDirAngle", refDirAngle);
  
  Environment::GetSingleton()->GetIntValue("VssTree.accessTimeThreshold", accessTimeThreshold);
  //= 1000;
  Environment::GetSingleton()->GetIntValue("VssTree.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("VssTree.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 VssTree split type "<<name<<endl;
		exit(1);
	  }
	
  Environment::GetSingleton()->GetBoolValue("VssTree.randomize", randomize);
  Environment::GetSingleton()->GetBoolValue("VssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
  Environment::GetSingleton()->GetBoolValue("VssTree.useRss", mUseRss);

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

  root = NULL;
  
  splitCandidates = new vector<SortableEntry>;
}


VssTree::~VssTree()
{
  if (root)
    delete root;
}




void
VssStatistics::Print(ostream &app) const
{
  app << "===== VssTree 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 VssTree statistics ==========\n";

}


void
VssTreeLeaf::UpdatePvsSize()
{
  if (!mValidPvs) {
	Intersectable::NewMail();
	int pvsSize = 0;
	for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
		ri != rays.end();
		ri++)
	  if ((*ri).mRay->IsActive()) {
		Intersectable *object;
#if BIDIRECTIONAL_RAY
		object = (*ri).mRay->mOriginObject;
		if (object && !object->Mailed()) {
		  pvsSize++;
		  object->Mail();
		}
#endif
		object = (*ri).mRay->mTerminationObject;
		if (object && !object->Mailed()) {
		  pvsSize++;
		  object->Mail();
		}
	  }
	mPvsSize = pvsSize;
	mValidPvs = true;

  }
}

bool
VssTree::ClipRay(
				 VssTreeNode::RayInfo &rayInfo,
				 const AxisAlignedBox3 &box
				 )
{
  float tmin, tmax;
  static Ray ray;
  ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
  box.ComputeMinMaxT(ray, &tmin, &tmax);
  if (tmin >= tmax)
	return false;
	
  if (tmin > 0.0f)
	rayInfo.SetMinT(tmin);

  if (tmax < 1.0f)
	rayInfo.SetMaxT(tmax);
	
  //  	vssRay->SetupEndPoints(
  //  												 origin,
  //  												 termination
  //  												 );
  return true;
}


void
VssTree::Construct(
				   VssRayContainer &rays,
				   AxisAlignedBox3 *forcedBoundingBox
				   )
{
  stat.Start();
  
  maxMemory = maxStaticMemory;

  if (root)
    delete root;
	
  root = new VssTreeLeaf(NULL, rays.size());
  // first construct a leaf that will get subdivide
  VssTreeLeaf *leaf = (VssTreeLeaf *) root;

  stat.nodes = 1;
  
  bbox.Initialize();
  dirBBox.Initialize();

  if (mUseRss)
	forcedBoundingBox = NULL;
	
  for(VssRayContainer::const_iterator ri = rays.begin();
      ri != rays.end();
      ri++) {
		
	VssTreeNode::RayInfo info(*ri);
	if (forcedBoundingBox)
	  if (!ClipRay(info, *forcedBoundingBox))
		continue;
	leaf->AddRay(info);
		
    bbox.Include((*ri)->GetOrigin());
    bbox.Include((*ri)->GetTermination());
    
		
	dirBBox.Include(Vector3(
							(*ri)->GetDirParametrization(0),
							(*ri)->GetDirParametrization(1),
							0
							)
					);
  }
	
	
  if ( forcedBoundingBox ) 
	bbox = *forcedBoundingBox;
	
  cout<<"Bbox = "<<bbox<<endl;
  cout<<"Dirr Bbox = "<<dirBBox<<endl;

  stat.rays = leaf->rays.size();
  leaf->UpdatePvsSize();
  leaf->ComputeEntropyImportance();

  stat.initialPvsSize = leaf->GetPvsSize();
  // Subdivide();
  root = Subdivide(TraversalData(leaf, bbox, 0));

  if (splitCandidates) {
    // force realease of this vector
    delete splitCandidates;
    splitCandidates = new vector<SortableEntry>;
  }
  
  stat.Stop();

  stat.Print(cout);
  cout<<"#Total memory="<<GetMemUsage()<<endl;

}

int
VssTree::UpdateSubdivision()
{
  priority_queue<TraversalData> tStack;
  //  stack<TraversalData> tStack;
  
  tStack.push(TraversalData(root, bbox, 0));
	
  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()) {
	  VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
										data.bbox,
										backBox,
										frontBox
										);
	  if (!node->IsLeaf()) {
		subdivided++;
		VssTreeInterior *interior = (VssTreeInterior *) 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 {
	  VssTreeInterior *interior = (VssTreeInterior *) 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;
}


VssTreeNode *
VssTree::Subdivide(const TraversalData &tdata)
{
  VssTreeNode *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();
    
    VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
									  data.bbox,
									  backBox,
									  frontBox
									  );
    if (result == NULL)
      result = node;
    
    if (!node->IsLeaf()) {
			
      VssTreeInterior *interior = (VssTreeInterior *) 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);
    }
  }

  return result;
}


// returns selected plane for subdivision
int
VssTree::SelectPlane(
					 VssTreeLeaf *leaf,
					 const AxisAlignedBox3 &box,
					 float &position,
					 int &raysBack,
					 int &raysFront,
					 int &pvsBack,
					 int &pvsFront
					 )
{

  int minDirDepth = 6;
  int axis;
  float costRatio;
	
  costRatio = BestCostRatio(leaf,
							axis,
							position,
							raysBack,
							raysFront,
							pvsBack,
							pvsFront
							);
#if DEBUG_SPLIT_COST
  cout<<axis<<" r="<<costRatio<<endl;
#endif
  if (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="<<axis<<endl;
#endif
	
  return axis;
}


float
VssTree::GetCostRatio(
					  VssTreeLeaf *leaf,
					  const int axis,
					  const float position,
					  const int raysBack,
					  const int raysFront,
					  const int pvsBack,
					  const int pvsFront
					  )
{
  bool costImportance = true;

  float ratio;
  AxisAlignedBox3 box;
  float minBox, maxBox;

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

  int pvsSize = leaf->GetPvsSize();

  if (!costImportance) {
	//		float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
	float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
	float newCost = ct_div_ci + sum/sizeBox;
	float oldCost = pvsSize;
	ratio = newCost/oldCost;
  } else {
	// importance based cost
#if 0
	float newContrib =
	  ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
	   (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
		
	//			float newContrib =
	//				sqr(pvsBack/(raysBack + Limits::Small)) +
	//				sqr(pvsFront/(raysFront + Limits::Small));
	float oldContrib = sqr(leaf->GetAvgRayContribution());
	ratio = oldContrib/newContrib;
#else
#if 1
	float newCost = raysBack*pvsBack  + raysFront*pvsFront;
	float oldCost = (float)leaf->rays.size()*pvsSize;
	ratio = newCost/oldCost;
#else
	float newCost = (pvsBack  + pvsFront)*0.5f;
	float oldCost = pvsSize;
	ratio = newCost/oldCost;
#endif
#endif
  }
	
  return ratio;
}
							

float
VssTree::EvalCostRatio(
					   VssTreeLeaf *leaf,
					   const int axis,
					   const float position,
					   int &raysBack,
					   int &raysFront,
					   int &pvsBack,
					   int &pvsFront
					   )
{
  raysBack = 0;
  raysFront = 0;
  pvsFront = 0;
  pvsBack = 0;

	
  Intersectable::NewMail(3);
	
  if (axis <= VssTreeNode::SPLIT_Z) {
    // this is the main ray classification loop!
    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		ri != leaf->rays.end();
		ri++)
      if ((*ri).mRay->IsActive()) {
		float t;
		// determine the side of this ray with respect to the plane
		int side = (*ri).ComputeRayIntersection(axis, position, t);
		//				(*ri).mRay->mSide = side;
				
		if (side <= 0)
		  raysBack++;
				
		if (side >= 0)
		  raysFront++;
				
		AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
      }
		
  } else {
		
	// directional split
    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		ri != leaf->rays.end();
		ri++)
      if ((*ri).mRay->IsActive()) {
				
		// determine the side of this ray with respect to the plane
		int side;
		if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
		  side = 1;
		else
		  side = -1;

		if (side <= 0)
		  raysBack++;
				
		if (side >= 0)
		  raysFront++;

		//				(*ri).mRay->mSide = side;
		AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
				
      }
  }

  float ratio = GetCostRatio(
							 leaf,
							 axis,
							 position,
							 raysBack,
							 raysFront,
							 pvsBack,
							 pvsFront);

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

  //	cout<<"ratio="<<ratio<<endl;

  return ratio;
}

float
VssTree::BestCostRatio(
					   VssTreeLeaf *leaf,
					   int &axis,
					   float &position,
					   int &raysBack,
					   int &raysFront,
					   int &pvsBack,
					   int &pvsFront
					   )
{
  int nRaysBack[6], nRaysFront[6];
  int nPvsBack[6], nPvsFront[6];
  float nPosition[6];
  float nCostRatio[6];
  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());
		
  
  for (axis = 0; axis < 5; axis++)
	if (mInterleaveDirSplits ||
		(axis < 3 && leaf->depth < mDirSplitDepth) ||
		(axis >= 3 && leaf->depth >= mDirSplitDepth)
		) {
	  if (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis) {
		
		
		if (splitType == ESplitRegular) {
		  if (axis < 3)
			nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
		  else
			nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
		  
		  nCostRatio[axis] = EvalCostRatio(leaf,
										   axis,
										   nPosition[axis],
										   nRaysBack[axis],
										   nRaysFront[axis],
										   nPvsBack[axis],
										   nPvsFront[axis]
										   );
		} else
		  if (splitType == ESplitHeuristic) {
			nCostRatio[axis] = EvalCostRatioHeuristic(
													  leaf,
													  axis,
													  nPosition[axis],
													  nRaysBack[axis],
													  nRaysFront[axis],
													  nPvsBack[axis],
													  nPvsFront[axis]);
		  } else
			if (splitType == ESplitHybrid) {
			  if (leaf->depth > 7)
				nCostRatio[axis] = EvalCostRatioHeuristic(
														  leaf,
														  axis,
														  nPosition[axis],
														  nRaysBack[axis],
														  nRaysFront[axis],
														  nPvsBack[axis],
														  nPvsFront[axis]);
			  else {
				if (axis < 3)
				  nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
				else
				  nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
				
				nCostRatio[axis] = EvalCostRatio(leaf,
												 axis,
												 nPosition[axis],
												 nRaysBack[axis],
												 nRaysFront[axis],
												 nPvsBack[axis],
												 nPvsFront[axis]
												 );
			  }
			} else {
			  cerr<<"VssTree: Unknown split heuristics\n";
			  exit(1);
			}
		
		
		if ( bestAxis == -1)
		  bestAxis = axis;
		else
		  if ( nCostRatio[axis] < nCostRatio[bestAxis] )
			bestAxis = axis;
	  }
	}
  
  axis = bestAxis;
  position = nPosition[bestAxis];

  raysBack = nRaysBack[bestAxis];
  raysFront = nRaysFront[bestAxis];

  pvsBack = nPvsBack[bestAxis];
  pvsFront = nPvsFront[bestAxis];
	
  return nCostRatio[bestAxis];
}

	
float
VssTree::EvalCostRatioHeuristic(
								VssTreeLeaf *leaf,
								const int axis,
								float &bestPosition,
								int &raysBack,
								int &raysFront,
								int &pvsBack,
								int &pvsFront
								)
{
  AxisAlignedBox3 box;
  float minBox, maxBox;
	
  if (axis < 3) {
	box = GetBBox(leaf);
	minBox = box.Min(axis);
	maxBox = box.Max(axis);
  } else {
	box = GetDirBBox(leaf);
	minBox = box.Min(axis-3);
	maxBox = box.Max(axis-3);
  }
	
  SortSubdivisionCandidates(leaf, axis);
  
  // go through the lists, count the number of objects left and right
  // and evaluate the following cost funcion:
  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
	
  int rl=0, rr = leaf->rays.size();
  int pl=0, pr = leaf->GetPvsSize();
  float sizeBox = maxBox - minBox;
	
  float minBand = minBox + 0.1*(maxBox - minBox);
  float maxBand = minBox + 0.9*(maxBox - minBox);
	
  float minRatio = 1e20;
	
  Intersectable::NewMail();
  // set all object as belonging to the fron pvs
  for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	  ri != leaf->rays.end();
	  ri++)
	if ((*ri).mRay->IsActive()) {
	  Intersectable *object = (*ri).mRay->mTerminationObject;
	  if (object)
		if (!object->Mailed()) {
		  object->Mail();
		  object->mCounter = 1;
		} else
		  object->mCounter++;
	}
	
  Intersectable::NewMail();
	
  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
	  ci < splitCandidates->end();
	  ci++) {
	VssRay *ray;
	switch ((*ci).type) {
	case SortableEntry::ERayMin: {
	  rl++;
	  ray = (VssRay *) (*ci).data;
	  Intersectable *object = ray->mTerminationObject;
	  if (object && !object->Mailed()) {
		object->Mail();
		pl++;
	  }
	  break;
	}
	case SortableEntry::ERayMax: {
	  rr--;
	  ray = (VssRay *) (*ci).data;
	  Intersectable *object = ray->mTerminationObject;
	  if (object) {
		if (--object->mCounter == 0)
		  pr--;
	  }
	  break;
	}
	}

	float position = (*ci).value;
		
	if (position > minBand && position < maxBand) {
			
	  float ratio = GetCostRatio(
								 leaf,
								 axis,
								 position,
								 rl,
								 rr,
								 pl,
								 pr);
			
			
	  //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	  //      cout<<"cost= "<<sum<<endl;
			
	  if (ratio < minRatio) {
		minRatio = ratio;
		bestPosition = position;
				
		raysBack = rl;
		raysFront = rr;
				
		pvsBack = pl;
		pvsFront = pr;
				
      }
    }
  }

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

void
VssTree::SortSubdivisionCandidates(
							 VssTreeLeaf *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(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
      ri < node->rays.end();
      ri++) {
	if ((*ri).mRay->IsActive()) {
	  if (axis < 3) {
		bool positive = (*ri).mRay->HasPosDir(axis);
		splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
												 SortableEntry::ERayMax,
												 (*ri).ExtrapOrigin(axis),
												 (void *)(*ri).mRay)
								   );
				
		splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
												 SortableEntry::ERayMin,
												 (*ri).ExtrapTermination(axis),
												 (void *)(*ri).mRay)
								   );
	  } else {
		float pos = (*ri).mRay->GetDirParametrization(axis-3);
		splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
												 pos - Limits::Small,
												 (void *)(*ri).mRay)
								   );
				
		splitCandidates->push_back(SortableEntry(SortableEntry::ERayMax,
												 pos + Limits::Small,
												 (void *)(*ri).mRay)
								   );
	  }
	}
  }

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


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

  // the node became a leaf -> evaluate stats for leafs
  VssTreeLeaf *leaf = (VssTreeLeaf *)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 (0 && leaf->GetAvgRayContribution() > termMaxRayContribution )
	stat.maxRayContribNodes++;
	
  if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
	stat.minSizeNodes++;
  }

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

}

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


VssTreeNode *
VssTree::SubdivideNode(
					   VssTreeLeaf *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;
  }
	
  float position;
	
  // first count ray sides
  int raysBack;
  int raysFront;
  int pvsBack;
  int pvsFront;
	
  // select subdivision axis
  int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
  //	cout<<axis<<" ";
	
  //	cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
  //		pvsFront<<endl;

  if (axis == -1) {
    return leaf;
  }
  
  stat.nodes+=2;
  stat.splits[axis]++;

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

  node->axis = axis;
  node->position = position;
  node->bbox = box;
  node->dirBBox = GetDirBBox(leaf);
  
  backBBox = box;
  frontBBox = box;

  VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
  VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);

  // replace a link from node's parent
  if (  leaf->parent )
    leaf->parent->ReplaceChildLink(leaf, node);
  // and setup child links
  node->SetupChildLinks(back, front);
	
  if (axis <= VssTreeNode::SPLIT_Z) {
	backBBox.SetMax(axis, position);
    frontBBox.SetMin(axis, position);
		
	for(VssTreeNode::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
		float t;
		int side = node->ComputeRayIntersection(*ri, t);
		
		if (side == 0) {
		  if ((*ri).mRay->HasPosDir(axis)) {
			back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
											  (*ri).mMinT,
											  t)
						 );
			front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
											   t,
											   (*ri).mMaxT));
		  } else {
			back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
											  t,
											  (*ri).mMaxT));
			front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
											   (*ri).mMinT,
											   t));
		  }
		} else
		  if (side == 1) 
			front->AddRay(*ri);
		  else
			back->AddRay(*ri);
      } else
		(*ri).mRay->Unref();
    }
  } else {
    // rays front/back
    
    
    for(VssTreeNode::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();

		int side;
		if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
		  side = 1;
		else
		  side = -1;
				
		if (side == 1)
		  front->AddRay(*ri);
		else
		  back->AddRay(*ri);
				
      } else
		(*ri).mRay->Unref();
    }
  }

  front->SetPvsSize(pvsFront);
  back->SetPvsSize(pvsBack);
  // compute entropy as well
  front->ComputeEntropyImportance();
  back->ComputeEntropyImportance();
  
  // update stats
  stat.rayRefs -= (int)leaf->rays.size();
  stat.rayRefs += raysBack + raysFront;

  
  delete leaf;
  return node;
}






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

  tstack.push(root);

  while (!tstack.empty()) {
    VssTreeNode *node = tstack.top();
    tstack.pop();
    
 
    if (!node->IsLeaf()) {
      VssTreeInterior *in = (VssTreeInterior *)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;
}




VssTreeNode *
VssTree::SubdivideLeaf(
					   VssTreeLeaf *leaf
					   )
{
  VssTreeNode *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
VssTree::UpdateRays(VssRayContainer &remove,
					VssRayContainer &add
					)
{
  VssTreeLeaf::NewMail();

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

  int inactive=0;

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


  //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
  
  for(VssRayContainer::const_iterator ri = add.begin();
      ri != add.end();
      ri++) {
	VssTreeNode::RayInfo info(*ri);
	if (ClipRay(info, bbox))
	  AddRay(info);
  }
}


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

  tstack.push(RayTraversalData(root, VssTreeNode::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...
      VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
      
      if (!leaf->Mailed()) {
		leaf->Mail();
		if (affectedLeaves)
		  affectedLeaves->push_back(leaf);
	
		if (removeAllScheduledRays) {
		  int tail = 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
VssTree::AddRay(VssTreeNode::RayInfo &info)
{

  stack<RayTraversalData> tstack;
  
  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...
      VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
      leaf->AddRay(data.rayData);
      stat.addedRayRefs++;
    }
  }
}

void
VssTree::TraverseInternalNode(
							  RayTraversalData &data,
							  stack<RayTraversalData> &tstack)
{
  VssTreeInterior *in =  (VssTreeInterior *) data.node;
  
  if (in->axis <= VssTreeNode::SPLIT_Z) {
	float t;
    // determine the side of this ray with respect to the plane
    int side = in->ComputeRayIntersection(data.rayData,
										  t);
    
  
    if (side == 0) {
      if (data.rayData.mRay->HasPosDir(in->axis)) {
		tstack.push(RayTraversalData(in->back,
									 VssTreeNode::RayInfo(data.rayData.mRay,
														  data.rayData.mMinT,
														  t))
					);
				
		tstack.push(RayTraversalData(in->front,
									 VssTreeNode::RayInfo(data.rayData.mRay,
														  t,
														  data.rayData.mMaxT
														  ))
					);
	
      } else {
		tstack.push(RayTraversalData(in->back,
									 VssTreeNode::RayInfo(data.rayData.mRay,
														  t,
														  data.rayData.mMaxT
														  ))
					);
				
		tstack.push(RayTraversalData(in->front,
									 VssTreeNode::RayInfo(data.rayData.mRay,
														  data.rayData.mMinT,
														  t))
					);
				
				
      }
    } else
      if (side == 1)
		tstack.push(RayTraversalData(in->front, data.rayData));
      else
		tstack.push(RayTraversalData(in->back, data.rayData));
  } 
  else {
    // directional split
	if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
	  tstack.push(RayTraversalData(in->front, data.rayData));
	else
	  tstack.push(RayTraversalData(in->back, data.rayData));
  }
}


int
VssTree::CollapseSubtree(VssTreeNode *sroot, const int time)
{
  // first count all rays in the subtree
  // use mail 1 for this purpose
  stack<VssTreeNode *> 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++;
    VssTreeNode *node = tstack.top();
    tstack.pop();

    if (node->IsLeaf()) {
      VssTreeLeaf *leaf = (VssTreeLeaf *) node;
      for(VssTreeNode::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(((VssTreeInterior *)node)->back);
      tstack.push(((VssTreeInterior *)node)->front);
    }
  }
  
  VssRay::NewMail();

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

  tstack.push( sroot );

  while (!tstack.empty()) {

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

    if (node->IsLeaf()) {
      VssTreeLeaf *leaf = (VssTreeLeaf *) node;
      
      for(VssTreeNode::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(((VssTreeInterior *)node)->back);
      tstack.push(((VssTreeInterior *)node)->front);
    }
  }
  
  // delete the node and all its children
  delete sroot;
  
  //   for(VssTreeNode::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.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
VssTree::GetPvsSize(const AxisAlignedBox3 &box) const
{
  stack<VssTreeNode *> tstack;
  tstack.push(root);

  Intersectable::NewMail();
  int pvsSize = 0;
	
  while (!tstack.empty()) {
    VssTreeNode *node = tstack.top();
    tstack.pop();
    
 
    if (node->IsLeaf()) {
	  VssTreeLeaf *leaf = (VssTreeLeaf *)node;
	  for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
		  ri != leaf->rays.end();
		  ri++)
		if ((*ri).mRay->IsActive()) {
		  Intersectable *object;
#if BIDIRECTIONAL_RAY
		  object = (*ri).mRay->mOriginObject;
		  if (object && !object->Mailed()) {
			pvsSize++;
			object->Mail();
		  }
#endif
		  object = (*ri).mRay->mTerminationObject;
		  if (object && !object->Mailed()) {
			pvsSize++;
			object->Mail();
		  }
		}
	} else {
	  VssTreeInterior *in = (VssTreeInterior *)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;
}

void
VssTree::GetRayContributionStatistics(
									  float &minRayContribution,
									  float &maxRayContribution,
									  float &avgRayContribution
									  )
{
  stack<VssTreeNode *> tstack;
  tstack.push(root);
	
  minRayContribution = 1.0f;
  maxRayContribution = 0.0f;
  float sumRayContribution = 0.0f;
  int leaves = 0;

  while (!tstack.empty()) {
    VssTreeNode *node = tstack.top();
    tstack.pop();
		
    if (node->IsLeaf()) {
	  leaves++;
	  VssTreeLeaf *leaf = (VssTreeLeaf *)node;
	  float c = leaf->GetImportance();
	  if (c > maxRayContribution)
		maxRayContribution = c;
	  if (c < minRayContribution)
		minRayContribution = c;
	  sumRayContribution += c;
			
	} else {
	  VssTreeInterior *in = (VssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }
	
  cout<<"sum="<<sumRayContribution<<endl;
  cout<<"leaves="<<leaves<<endl;
  avgRayContribution = sumRayContribution/(float)leaves;
}


int
VssTree::GenerateRays(const float ratioPerLeaf,
					  SimpleRayContainer &rays)
{
  stack<VssTreeNode *> tstack;
  tstack.push(root);
	
  while (!tstack.empty()) {
    VssTreeNode *node = tstack.top();
    tstack.pop();
		
    if (node->IsLeaf()) {
	  VssTreeLeaf *leaf = (VssTreeLeaf *)node;
	  float c = leaf->GetImportance();
	  int num = (c*ratioPerLeaf + 0.5);
	  //			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));
	  }
			
	} else {
	  VssTreeInterior *in = (VssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }

  return rays.size();
}

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

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


void
VssTree::GenerateLeafRays(VssTreeLeaf *leaf,
						  const int numberOfRays,
						  SimpleRayContainer &rays)
{
  int nrays = (int)leaf->rays.size();
  for (int i=0; i < numberOfRays; i++) {
	// pickup 3 random rays
	int r1 = (int)RandomValue(0, nrays-1);
	int r2 = (int)RandomValue(0, nrays-1);
	int r3 = (int)RandomValue(0, nrays-1);
		
	Vector3 o1 = leaf->rays[r1].Extrap(RandomValue(leaf->rays[r1].GetMinT(), 
												   leaf->rays[r1].GetMaxT()));

	Vector3 o2 = leaf->rays[r2].Extrap(RandomValue(leaf->rays[r2].GetMinT(), 
												   leaf->rays[r2].GetMaxT()));
		
	Vector3 o3 = leaf->rays[r3].Extrap(RandomValue(leaf->rays[r3].GetMinT(), 
												   leaf->rays[r3].GetMaxT()));

	const float overlap = 0.1f;
		
	Vector3 origin, direction;
	bool useExtendedConvexCombination = true;
	if (useExtendedConvexCombination) {
	  float w1, w2, w3;
	  GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
	  origin = w1*o1 + w2*o2 + w3*o3;
	  direction =
		w1*leaf->rays[r1].mRay->GetDir() +
		w2*leaf->rays[r2].mRay->GetDir() +
		w3*leaf->rays[r3].mRay->GetDir();
	} else {
	  origin = GetBBox(leaf).GetRandomPoint();
	  Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
	  direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
	}
	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	rays.push_back(SimpleRay(origin, direction));
  }
}

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

  vector<VssTreeLeaf *> leaves;
	
  CollectLeaves(leaves);
	
  sort(leaves.begin(),
	   leaves.end(),
	   GreaterContribution);
			 

  float sumContrib = 0.0;
  int i;
  int k = 0;
  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	if (ValidLeaf(leaves[i])) {
	  float c = leaves[i]->GetImportance();
	  sumContrib += c;
	  //		cout<<"ray contrib "<<i<<" : "<<c<<endl;
	  k++;
	}
				
  float avgContrib = sumContrib/numberOfLeaves;
  float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
  k = 0;
  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	if (ValidLeaf(leaves[i])) {
	  k++;
	  VssTreeLeaf *leaf = leaves[i];
	  float c = leaf->GetImportance();
	  int num = (int)(c*ratioPerLeaf + 0.5f);
	  GenerateLeafRays(leaf, num, rays);
	}

  return (int)rays.size();
}


float
VssTree::GetAvgPvsSize()
{
  stack<VssTreeNode *> tstack;
  tstack.push(root);

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


  return sumPvs/(float)leaves;
}

	
float
VssTreeLeaf::GetImportance()  const
{

  if (1) {
	//	return GetAvgRayContribution();
	return (float)GetPvsSize();
  } else {
	// return GetAvgRayContribution()*mEntropyImportance;
	//return GetAvgRayContribution();
	return mEntropyImportance;
  }
}


float
VssTreeLeaf::ComputePvsEntropy()
{
  int samples = 0;
  Intersectable::NewMail();
  // set all object as belonging to the fron pvs
  for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
	  ri != rays.end();
	  ri++)
	if ((*ri).mRay->IsActive()) {
	  Intersectable *object = (*ri).mRay->mTerminationObject;
	  if (object) {
		if (!object->Mailed()) {
		  object->Mail();
		  object->mCounter = 1;
		} else
		  object->mCounter++;
		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).mRay->mTerminationObject;
		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
VssTreeLeaf::ComputeRayLengthEntropy()
{
  // get sum of all ray lengths
  // consider only passing rays or originating rays
  float sum = 0.0f;
  int samples = 0;

  for(RayInfoContainer::const_iterator ri = rays.begin();
      ri != rays.end();
      ri++) {
	int rayClass = (*ri).GetRayClass();
	if (
		rayClass == RayInfo::PASSING_RAY
		//		||
		//		rayClass == RayInfo::SOURCE_RAY
		//		rayClass == RayInfo::TERMINATION_RAY
		) {
	  float c = 1.0f - (*ri).GetMaxT();
	  sum += (*ri).mRay->GetSize()*c;
	  samples++;
	}
  }
  
  float entropy = 0.0f;
  
  if (samples > 1) {
	for(RayInfoContainer::const_iterator ri = rays.begin();
		ri != rays.end();
		ri++) {
	  int rayClass = (*ri).GetRayClass();
	  if (
		  rayClass == RayInfo::PASSING_RAY
		  //		  ||
		  //		  rayClass == RayInfo::SOURCE_RAY
		  // rayClass == RayInfo::TERMINATION_RAY
		  ) {
		float c = 1.0f - (*ri).GetMaxT();
		float p = (*ri).mRay->GetSize()*c/sum;
		entropy -= p*log(p);
	  }
	}
	entropy = entropy/log((float)samples);
  } else
	entropy = 1.0f;
  
  return entropy;
}
  
float
VssTreeLeaf::ComputeRayTerminationEntropy()
{

  return 0.0f;
}


void
VssTreeLeaf::ComputeEntropyImportance()
{
  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
  mEntropyImportance = 1.0f - ComputePvsEntropy();
  
  //  cout<<"ei="<<mEntropyImportance<<" ";
}


int
VssTree::CollectRays(VssRayContainer &rays,
					 const int number)
{
  VssRayContainer allRays;
  CollectRays(allRays);
  
  int desired = min(number, (int)allRays.size());
  float prob = desired/(float)allRays.size();
  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 rays.size();
}


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

  stack<VssTreeNode *> tstack;
  tstack.push(root);
  
  while (!tstack.empty()) {
    VssTreeNode *node = tstack.top();
    tstack.pop();
	
    if (node->IsLeaf()) {
	  VssTreeLeaf *leaf = (VssTreeLeaf *)node;
	  // update pvs size
	  VssTreeNode::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 {
	  VssTreeInterior *in = (VssTreeInterior *)node;
	  // both nodes for directional splits
	  tstack.push(in->front);
	  tstack.push(in->back);
	}
  }
  
  return rays.size();
}

}