Index: trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj	(revision 406)
+++ trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj	(revision 408)
@@ -291,8 +291,8 @@
 			</File>
 			<File
-				RelativePath="..\src\ViewCellKd.cpp">
-			</File>
-			<File
-				RelativePath="..\src\ViewCellKd.h">
+				RelativePath="..\src\VspKdTree.cpp">
+			</File>
+			<File
+				RelativePath="..\src\VspKdTree.h">
 			</File>
 			<File
Index: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp	(revision 408)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp	(revision 408)
@@ -0,0 +1,1492 @@
+
+// ================================================================
+// $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 "VspKdTree.h"
+
+#include "Environment.h"
+#include "VssRay.h"
+#include "Intersectable.h"
+#include "Ray.h"
+
+
+#define DEBUG_DIR_SPLIT 0
+
+
+
+// Static variables
+int
+VspKdTreeLeaf::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
+VspKdTree::VspKdTree()
+{
+  environment->GetIntValue("VspKdTree.maxDepth", termMaxDepth);
+	environment->GetIntValue("VspKdTree.minPvs", termMinPvs);
+	environment->GetIntValue("VspKdTree.minRays", termMinRays);
+	environment->GetFloatValue("VspKdTree.maxRayContribution", termMaxRayContribution);
+  environment->GetFloatValue("VspKdTree.maxCostRatio", termMaxCostRatio);
+
+  environment->GetFloatValue("VspKdTree.minSize", termMinSize);
+  termMinSize = sqr(termMinSize);
+	
+  environment->GetFloatValue("VspKdTree.refDirBoxMaxSize", refDirBoxMaxSize);
+  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
+  
+  environment->GetFloatValue("VspKdTree.epsilon", epsilon);
+  environment->GetFloatValue("VspKdTree.ct_div_ci", ct_div_ci);
+	
+  environment->GetFloatValue("VspKdTree.maxTotalMemory", maxTotalMemory);
+  environment->GetFloatValue("VspKdTree.maxStaticMemory", maxStaticMemory);
+  
+  
+
+
+  float refDirAngle;
+  environment->GetFloatValue("VspKdTree.refDirAngle", refDirAngle);
+
+  environment->GetIntValue("VspKdTree.accessTimeThreshold", accessTimeThreshold);
+  //= 1000;
+  environment->GetIntValue("VspKdTree.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->GetStringValue("VspKdTree.splitType", sname);
+  string name(sname);
+	
+  if (name.compare("regular") == 0)
+    splitType = ESplitRegular;
+  else
+    if (name.compare("heuristic") == 0)
+      splitType = ESplitHeuristic;
+		else {
+			cerr<<"Invalid VspKdTree split type "<<name<<endl;
+			exit(1);
+		}
+	
+  environment->GetBoolValue("VspKdTree.randomize", randomize);
+
+  root = NULL;
+  
+  splitCandidates = new vector<SortableEntry>;
+}
+
+
+VspKdTree::~VspKdTree()
+{
+  if (root)
+    delete root;
+}
+
+
+
+
+void
+VspKdStatistics::Print(ostream &app) const
+{
+  app << "===== VspKdTree 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 minCost )\n"<<
+    minPvsNodes*100/(double)Leaves()<<endl;
+
+  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\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_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 VspKdTree statistics ==========\n";
+
+}
+
+
+void
+VspKdTreeLeaf::UpdatePvsSize()
+{
+	if (!mValidPvs) {
+		Intersectable::NewMail();
+		int pvsSize = 0;
+		for(VspKdTreeNode::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;
+	}
+}
+
+
+void
+VspKdTree::Construct(
+									 VssRayContainer &rays,
+									 AxisAlignedBox3 *forcedBoundingBox
+									 )
+{
+  stat.Start();
+  
+	maxMemory = maxStaticMemory;
+
+  if (root)
+    delete root;
+
+  root = new VspKdTreeLeaf(NULL, rays.size());
+  // first construct a leaf that will get subdivide
+  VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) root;
+
+  stat.nodes = 1;
+  
+  bbox.Initialize();
+  dirBBox.Initialize();
+  
+  for(VssRayContainer::const_iterator ri = rays.begin();
+      ri != rays.end();
+      ri++) {
+    leaf->AddRay(VspKdTreeNode::RayInfo(*ri));
+
+    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();
+  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;
+
+}
+
+
+
+VspKdTreeNode *
+VspKdTree::Subdivide(const TraversalData &tdata)
+{
+  VspKdTreeNode *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();
+    
+    VspKdTreeNode *node = SubdivideNode((VspKdTreeLeaf *) data.node,
+																			data.bbox,
+																			backBox,
+																			frontBox
+																			);
+    if (result == NULL)
+      result = node;
+    
+    if (!node->IsLeaf()) {
+			
+      VspKdTreeInterior *interior = (VspKdTreeInterior *) 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
+VspKdTree::SelectPlane(
+										 VspKdTreeLeaf *leaf,
+										 const AxisAlignedBox3 &box,
+										 float &position,
+										 int &raysBack,
+										 int &raysFront,
+										 int &pvsBack,
+										 int &pvsFront
+										 )
+{
+
+  int minDirDepth = 6;
+  int axis;
+	float costRatio;
+	
+  if (splitType == ESplitRegular) {
+		costRatio = BestCostRatioRegular(leaf,
+																		 axis,
+																		 position,
+																		 raysBack,
+																		 raysFront,
+																		 pvsBack,
+																		 pvsFront
+																		 );
+
+	} else {
+    if (splitType == ESplitHeuristic)
+      costRatio = BestCostRatioHeuristic(leaf,
+																				 axis,
+																				 position,
+																				 raysBack,
+																				 raysFront,
+																				 pvsBack,
+																				 pvsFront
+																				 );
+		else {
+			cerr<<"VspKdTree: Unknown split heuristics\n";
+			exit(1);
+		}
+  }
+
+	if (costRatio > termMaxCostRatio) {
+		//		cout<<"Too big cost ratio "<<costRatio<<endl;
+		return -1;
+	}
+
+#if 0	
+	cout<<
+		"pvs="<<leaf->mPvsSize<<
+		" rays="<<leaf->rays.size()<<
+		" rc="<<leaf->GetAvgRayContribution()<<
+		" axis="<<axis<<endl;
+#endif
+	
+	return axis;
+}
+
+
+							
+
+float
+VspKdTree::EvalCostRatio(
+											 VspKdTreeLeaf *leaf,
+											 const int axis,
+											 const float position,
+											 int &raysBack,
+											 int &raysFront,
+											 int &pvsBack,
+											 int &pvsFront
+											 )
+{
+	raysBack = 0;
+	raysFront = 0;
+	pvsFront = 0;
+	pvsBack = 0;
+
+	float newCost;
+
+	Intersectable::NewMail(3);
+	
+	// eval pvs size
+	int pvsSize = leaf->GetPvsSize();
+
+	Intersectable::NewMail(3);
+
+	if (axis <= VspKdTreeNode::SPLIT_Z) {
+    // this is the main ray classification loop!
+    for(VspKdTreeNode::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 = (*ri).ComputeRayIntersection(axis, position, (*ri).mRay->mT);
+				//				(*ri).mRay->mSide = side;
+				
+				if (side <= 0)
+					raysBack++;
+				
+				if (side >= 0)
+					raysFront++;
+				
+				AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+      }
+		
+		AxisAlignedBox3 box = GetBBox(leaf);
+		
+		float minBox = box.Min(axis);
+		float maxBox = box.Max(axis);
+		float sizeBox = maxBox - minBox;
+		
+		//		float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+		float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
+
+		newCost = ct_div_ci + sum/sizeBox;
+		
+  } else {
+		
+		// directional split
+    for(VspKdTreeNode::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);
+
+      }
+		
+		AxisAlignedBox3 box = GetDirBBox(leaf);
+		
+		float minBox = box.Min(axis-3);
+		float maxBox = box.Max(axis-3);
+		float sizeBox = maxBox - minBox;
+		
+		//		float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+		float sum =
+			pvsBack*(position - minBox) +
+			pvsFront*(maxBox - position);
+		//		float sum = pvsBack + pvsFront;
+		newCost = ct_div_ci + sum/sizeBox;
+  }
+
+	//	cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
+	//  float oldCost = leaf->rays.size();
+	float oldCost = pvsSize;
+  float ratio = newCost/oldCost;
+	//	cout<<"ratio="<<ratio<<endl;
+	return ratio;
+}
+
+float
+VspKdTree::BestCostRatioRegular(
+															VspKdTreeLeaf *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;
+
+	bool onlyDrivingAxis = true; 
+
+	for (axis = 0; axis < 6; axis++) {
+		if (!onlyDrivingAxis || axis == sAxis || axis == dAxis) {
+			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]
+																			 );
+			
+			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
+VspKdTree::BestCostRatioHeuristic(
+																VspKdTreeLeaf *leaf,
+																int &axis,
+																float &position,
+																int &raysBack,
+																int &raysFront,
+																int &pvsBack,
+																int &pvsFront
+																)
+{
+	AxisAlignedBox3 box = GetBBox(leaf);
+	//	AxisAlignedBox3 dirBox = GetDirBBox(node);
+	
+	axis = box.Size().DrivingAxis();
+		
+	SortSplitCandidates(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 minBox = box.Min(axis);
+	float maxBox = box.Max(axis);
+	float sizeBox = maxBox - minBox;
+	
+	float minBand = minBox + 0.1*(maxBox - minBox);
+	float maxBand = minBox + 0.9*(maxBox - minBox);
+	
+	float sum = rr*sizeBox;
+	float minSum = 1e20;
+	
+	Intersectable::NewMail();
+	// set all object as belonging to the fron pvs
+	for(VspKdTreeNode::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;
+		}
+		}
+		if ((*ci).value > minBand && (*ci).value < maxBand) {
+			
+			sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);
+			
+			//      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
+			//      cout<<"cost= "<<sum<<endl;
+			
+			if (sum < minSum) {
+				minSum = sum;
+				position = (*ci).value;
+				
+				raysBack = rl;
+				raysFront = rr;
+				
+				pvsBack = pl;
+				pvsFront = pr;
+				
+      }
+    }
+  }
+  
+  float oldCost = leaf->GetPvsSize();
+  float newCost = ct_div_ci + minSum/sizeBox;
+  float ratio = newCost/oldCost;
+  
+  //  cout<<"===================="<<endl;
+  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
+  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
+	return ratio;
+}
+
+void
+VspKdTree::SortSplitCandidates(
+														 VspKdTreeLeaf *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(VspKdTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
+      ri < node->rays.end();
+      ri++) {
+    bool positive = (*ri).mRay->HasPosDir(axis);
+    splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
+																						 SortableEntry::ERayMax,
+																						 (*ri).ExtrapOrigin(axis),
+																						 (void *)&*ri)
+															 );
+		
+    splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
+																						 SortableEntry::ERayMin,
+																						 (*ri).ExtrapTermination(axis),
+																						 (void *)&*ri)
+															 );
+  }
+	
+  stable_sort(splitCandidates->begin(), splitCandidates->end());
+}
+
+
+void
+VspKdTree::EvaluateLeafStats(const TraversalData &data)
+{
+
+  // the node became a leaf -> evaluate stats for leafs
+  VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)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 = leaf->rays.size();
+
+}
+
+
+
+VspKdTreeNode *
+VspKdTree::SubdivideNode(
+											 VspKdTreeLeaf *leaf,
+											 const AxisAlignedBox3 &box,
+											 AxisAlignedBox3 &backBBox,
+											 AxisAlignedBox3 &frontBBox
+											 )
+{
+  
+  if ( (leaf->GetPvsSize() < termMinPvs) ||
+			 (leaf->rays.size() < termMinRays) ||
+			 //			 (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
+       (leaf->depth >= termMaxDepth) ||
+			 SqrMagnitude(box.Size()) <= termMinSize ) {
+
+#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<<"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
+  VspKdTreeInterior *node = new VspKdTreeInterior(leaf->parent);
+
+  node->axis = axis;
+  node->position = position;
+  node->bbox = box;
+  node->dirBBox = GetDirBBox(leaf);
+  
+  backBBox = box;
+  frontBBox = box;
+
+  VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
+	back->SetPvsSize(pvsBack);
+  VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
+	front->SetPvsSize(pvsFront);
+
+  // 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 <= VspKdTreeNode::SPLIT_Z) {
+		backBBox.SetMax(axis, position);
+    frontBBox.SetMin(axis, position);
+		
+		for(VspKdTreeNode::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();
+
+				// determine the side of this ray with respect to the plane
+				int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
+
+				if (side == 0) {
+					if ((*ri).mRay->HasPosDir(axis)) {
+						back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+																							(*ri).mMinT,
+																							(*ri).mRay->mT)
+												 );
+						front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+																							 (*ri).mRay->mT,
+																							 (*ri).mMaxT));
+					} else {
+						back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+																							(*ri).mRay->mT,
+																							(*ri).mMaxT));
+						front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+																							 (*ri).mMinT,
+																							 (*ri).mRay->mT));
+					}
+				} else
+					if (side == 1) 
+						front->AddRay(*ri);
+					else
+						back->AddRay(*ri);
+      } else
+				(*ri).mRay->Unref();
+    }
+  } else {
+    // rays front/back
+    
+    
+    for(VspKdTreeNode::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();
+    }
+  }
+	
+  // update stats
+  stat.rayRefs -= leaf->rays.size();
+  stat.rayRefs += raysBack + raysFront;
+
+  
+  delete leaf;
+  return node;
+}
+
+
+
+
+
+
+int
+VspKdTree::ReleaseMemory(const int time)
+{
+  stack<VspKdTreeNode *> tstack;
+  
+  // find a node in the tree which subtree will be collapsed
+  int maxAccessTime = time - accessTimeThreshold;
+  int released;
+
+  tstack.push(root);
+
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
+    if (!node->IsLeaf()) {
+      VspKdTreeInterior *in = (VspKdTreeInterior *)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;
+}
+
+
+
+
+VspKdTreeNode *
+VspKdTree::SubdivideLeaf(
+											 VspKdTreeLeaf *leaf,
+											 const float sizeThreshold
+											 )
+{
+  VspKdTreeNode *node = leaf;
+
+  AxisAlignedBox3 leafBBox = GetBBox(leaf);
+
+	static int pass = 0;
+	pass ++;
+	
+  // check if we should perform a dynamic subdivision of the leaf
+  if (
+			//      leaf->rays.size() > (unsigned)termMinCost &&
+			(leaf->GetPvsSize() >= termMinPvs) &&
+      SqrMagnitude(leafBBox.Size()) > sizeThreshold) {
+    
+		// memory check and realese...
+    if (GetMemUsage() > maxTotalMemory) {
+      ReleaseMemory( pass );
+    }
+    
+    AxisAlignedBox3 backBBox, frontBBox;
+
+    // subdivide the node
+    node =
+      SubdivideNode(leaf,
+										leafBBox,
+										backBBox,
+										frontBBox
+										);
+  }
+  return node;
+}
+
+
+
+void
+VspKdTree::UpdateRays(VssRayContainer &remove,
+										VssRayContainer &add
+										)
+{
+  VspKdTreeLeaf::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++) {
+    AddRay(*ri);
+  }
+  
+
+}
+
+
+void
+VspKdTree::RemoveRay(VssRay *ray,
+									 vector<VspKdTreeLeaf *> *affectedLeaves,
+									 const bool removeAllScheduledRays
+									 )
+{
+	
+  stack<RayTraversalData> tstack;
+
+  tstack.push(RayTraversalData(root, VspKdTreeNode::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...
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)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
+VspKdTree::AddRay(VssRay *ray)
+{
+
+  stack<RayTraversalData> tstack;
+  
+  tstack.push(RayTraversalData(root, VspKdTreeNode::RayInfo(ray)));
+  
+  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...
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)data.node;
+      leaf->AddRay(data.rayData);
+      stat.addedRayRefs++;
+    }
+  }
+}
+
+void
+VspKdTree::TraverseInternalNode(
+															RayTraversalData &data,
+															stack<RayTraversalData> &tstack)
+{
+  VspKdTreeInterior *in =  (VspKdTreeInterior *) data.node;
+  
+  if (in->axis <= VspKdTreeNode::SPLIT_Z) {
+
+    // determine the side of this ray with respect to the plane
+    int side = in->ComputeRayIntersection(data.rayData,
+																					data.rayData.mRay->mT);
+    
+  
+    if (side == 0) {
+      if (data.rayData.mRay->HasPosDir(in->axis)) {
+				tstack.push(RayTraversalData(in->back,
+																		 VspKdTreeNode::RayInfo(data.rayData.mRay,
+																													data.rayData.mMinT,
+																													data.rayData.mRay->mT))
+										);
+				
+				tstack.push(RayTraversalData(in->front,
+																		 VspKdTreeNode::RayInfo(data.rayData.mRay,
+																													data.rayData.mRay->mT,
+																													data.rayData.mMaxT
+																													))
+										);
+	
+      } else {
+				tstack.push(RayTraversalData(in->back,
+																		 VspKdTreeNode::RayInfo(data.rayData.mRay,
+																													data.rayData.mRay->mT,
+																													data.rayData.mMaxT
+																													))
+										);
+				
+				tstack.push(RayTraversalData(in->front,
+																		 VspKdTreeNode::RayInfo(data.rayData.mRay,
+																													data.rayData.mMinT,
+																													data.rayData.mRay->mT))
+										);
+				
+				
+      }
+    } 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
+VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
+{
+  // first count all rays in the subtree
+  // use mail 1 for this purpose
+  stack<VspKdTreeNode *> 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++;
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+
+    if (node->IsLeaf()) {
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
+      for(VspKdTreeNode::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(((VspKdTreeInterior *)node)->back);
+      tstack.push(((VspKdTreeInterior *)node)->front);
+    }
+  }
+  
+  VssRay::NewMail();
+
+  // create a new node that will hold the rays
+  VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf( sroot->parent, rayCount );
+  if (  newLeaf->parent )
+    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
+  
+
+  tstack.push( sroot );
+
+  while (!tstack.empty()) {
+
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+
+    if (node->IsLeaf()) {
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
+      
+      for(VspKdTreeNode::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(((VspKdTreeInterior *)node)->back);
+      tstack.push(((VspKdTreeInterior *)node)->front);
+    }
+  }
+  
+  // delete the node and all its children
+  delete sroot;
+  
+//   for(VspKdTreeNode::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
+VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
+{
+	stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+
+	Intersectable::NewMail();
+	int pvsSize = 0;
+	
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
+    if (node->IsLeaf()) {
+			VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+			for(VspKdTreeNode::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 {
+			VspKdTreeInterior *in = (VspKdTreeInterior *)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
+VspKdTree::GetRayContributionStatistics(
+																			float &minRayContribution,
+																			float &maxRayContribution,
+																			float &avgRayContribution
+																			)
+{
+	stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+	
+	minRayContribution = 1.0f;
+	maxRayContribution = 0.0f;
+	float sumRayContribution = 0.0f;
+	int leaves = 0;
+
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+			leaves++;
+			VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+			float c = leaf->GetAvgRayContribution();
+			if (c > maxRayContribution)
+				maxRayContribution = c;
+			if (c < minRayContribution)
+				minRayContribution = c;
+			sumRayContribution += c;
+			
+		} else {
+			VspKdTreeInterior *in = (VspKdTreeInterior *)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
+VspKdTree::GenerateRays(const float ratioPerLeaf,
+											SimpleRayContainer &rays)
+{
+	stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+	
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+			VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+			float c = leaf->GetAvgRayContribution();
+			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 = 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));
+			}
+			
+		} else {
+			VspKdTreeInterior *in = (VspKdTreeInterior *)node;
+			// both nodes for directional splits
+			tstack.push(in->front);
+			tstack.push(in->back);
+		}
+	}
+
+	return rays.size();
+}
+
+
+float
+VspKdTree::GetAvgPvsSize()
+{
+	stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+
+	int sumPvs = 0;
+	int leaves = 0;
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+			VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+			// update pvs size
+			leaf->UpdatePvsSize();
+			sumPvs += leaf->GetPvsSize();
+			leaves++;
+		} else {
+			VspKdTreeInterior *in = (VspKdTreeInterior *)node;
+			// both nodes for directional splits
+			tstack.push(in->front);
+			tstack.push(in->back);
+		}
+	}
+
+
+	return sumPvs/(float)leaves;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 406)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 408)
@@ -1322,5 +1322,40 @@
 
 
-	RegisterOption("VssTree.maxDepth", optInt, "kd_depth=", "12");
+	RegisterOption("VspKdTree.maxDepth", optInt, "kd_depth=", "12");
+  RegisterOption("VspKdTree.minPvs", optInt, "kd_minpvs=", "1");
+  RegisterOption("VspKdTree.minRays", optInt, "kd_minrays=", "10");
+	RegisterOption("VspKdTree.maxCostRatio", optFloat, "maxcost=", "0.95");
+	RegisterOption("VspKdTree.maxRayContribution", optFloat, "maxraycontrib=", "0.5");
+
+	RegisterOption("VspKdTree.epsilon", optFloat, "kd_eps=", "1e-6");
+  RegisterOption("VspKdTree.ct_div_ci", optFloat, "kd_ctdivci=", "1.0");
+  RegisterOption("VspKdTree.randomize", optBool, "randomize", "false");
+  RegisterOption("VspKdTree.splitType", optString, "split=", "queries");
+  RegisterOption("VspKdTree.numberOfEndPointDomains", optInt, "endpoints=", "10000");
+
+  RegisterOption("VspKdTree.minSize", optFloat, "minsize=", "0.001");
+
+  RegisterOption("VspKdTree.maxTotalMemory", optFloat, "mem=", "60.0");
+  RegisterOption("VspKdTree.maxStaticMemory", optFloat, "statmem=", "8.0");
+
+  RegisterOption("VspKdTree.queryType", optString, "qtype=", "static");
+
+  RegisterOption("VspKdTree.queryPosWeight", optFloat, "qposweight=", "0.0");
+  RegisterOption("VspKdTree.useRefDirSplits", optBool, "refdir", "false");
+  RegisterOption("VspKdTree.refDirAngle", optFloat, "refangle=", "10");
+  RegisterOption("VspKdTree.refDirBoxMaxSize", optFloat, "refboxsize=", "0.1");
+  RegisterOption("VspKdTree.accessTimeThreshold", optInt, "accesstime=", "1000");
+  RegisterOption("VspKdTree.minCollapseDepth", optInt, "colldepth=", "4");
+
+
+
+
+
+
+
+
+
+
+  	RegisterOption("VssTree.maxDepth", optInt, "kd_depth=", "12");
   RegisterOption("VssTree.minPvs", optInt, "kd_minpvs=", "1");
   RegisterOption("VssTree.minRays", optInt, "kd_minrays=", "10");
@@ -1347,5 +1382,4 @@
   RegisterOption("VssTree.accessTimeThreshold", optInt, "accesstime=", "1000");
   RegisterOption("VssTree.minCollapseDepth", optInt, "colldepth=", "4");
-
 }
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h	(revision 408)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h	(revision 408)
@@ -0,0 +1,811 @@
+// ================================================================
+// $Id$
+// ****************************************************************
+// 
+// Initial coding by
+/**
+   @author Jiri Bittner
+*/
+
+#ifndef __VSPKDTREE_H__
+#define __VSPKDTREE_H__
+
+// Standard headers
+#include <iomanip>
+#include <vector>
+#include <functional>
+#include <stack>
+
+
+// User headers
+#include "VssRay.h"
+#include "AxisAlignedBox3.h"
+
+
+#define USE_KDNODE_VECTORS 1
+#define _RSS_STATISTICS
+#define _RSS_TRAVERSAL_STATISTICS
+
+
+#include "Statistics.h"
+#include "Ray.h"
+
+// --------------------------------------------------------------
+// Static statistics for kd-tree search
+// --------------------------------------------------------------
+class VspKdStatistics:  public StatisticsBase
+{
+public:  
+  // total number of nodes
+  int nodes;
+  // number of splits along each of the axes
+  int splits[7];
+  // totals number of rays
+  int rays;
+	// initial size of the pvs
+  int initialPvsSize;
+	// total number of query domains
+  int queryDomains;
+  // total number of ray references
+  int rayRefs;
+
+	// max depth nodes
+  int maxDepthNodes;
+  // max depth nodes
+  int minPvsNodes;
+  int minRaysNodes;
+	// max ray contribution nodes
+  int maxRayContribNodes;
+	// max depth nodes
+  int minSizeNodes;
+	
+  // max number of rays per node
+  int maxRayRefs;
+  // number of dynamically added ray refs
+  int addedRayRefs;
+  // number of dynamically removed ray refs
+  int removedRayRefs;
+  
+  // Constructor
+  VspKdStatistics() {
+    Reset();
+  }
+
+  int Nodes() const {return nodes;}
+  int Interior() const { return nodes/2; }
+  int Leaves() const { return (nodes/2) + 1; }
+
+  void Reset() {
+    nodes = 0;
+    for (int i=0; i<7; i++)
+      splits[i] = 0;
+    rays = queryDomains = 0;
+    rayRefs = 0;
+    maxDepthNodes = 0;
+    minPvsNodes = 0;
+    minRaysNodes = 0;
+    maxRayRefs = 0;
+    addedRayRefs = removedRayRefs = 0;
+		initialPvsSize = 0;
+		maxRayContribNodes = 0;
+		minSizeNodes = 0;
+  }
+
+  
+  void
+  Print(ostream &app) const;
+
+  friend ostream &operator<<(ostream &s, const VspKdStatistics &stat) {
+    stat.Print(s);
+    return s;
+  }
+  
+};
+
+
+// --------------------------------------------------------------
+// For sorting rays
+// --------------------------------------------------------------
+struct  SortableEntry
+{
+  enum EType {
+    ERayMin,
+    ERayMax
+  };
+
+  int type;
+  float value;
+  void *data;
+  
+  SortableEntry() {}
+  SortableEntry(const int t, const float v, void *d):type(t),
+																										 value(v),
+																										 data(d) {}
+	
+  friend bool operator<(const SortableEntry &a, const SortableEntry &b) {
+    return a.value < b.value;
+  }
+};
+
+
+class VspKdTreeInterior;
+
+
+// --------------------------------------------------------------
+// KD-tree node - abstract class
+// --------------------------------------------------------------
+class VspKdTreeNode
+{
+public:
+
+#define USE_FIXEDPOINT_T 1
+
+struct RayInfo
+{
+	// pointer to the actual ray
+	VssRay *mRay;
+	
+	// endpoints  - do we need them?
+#if USE_FIXEDPOINT_T
+	short mMinT, mMaxT;
+#else
+	float mMinT, mMaxT;
+#endif
+	
+	RayInfo():mRay(NULL) {}
+	
+	RayInfo(VssRay *r):mRay(r), mMinT(0),
+#if USE_FIXEDPOINT_T
+#define FIXEDPOINT_ONE 0x7FFE
+												//			  mMaxT(0xFFFF)
+																			mMaxT(FIXEDPOINT_ONE)
+#else
+		mMaxT(1.0f)
+#endif
+	{}
+	
+	RayInfo(VssRay *r,
+										const float _min,
+										const float _max
+										):mRay(r) {
+		SetMinT(_min);
+		SetMaxT(_max);
+	}
+	
+	RayInfo(VssRay *r,
+										const short _min,
+										const float _max
+										):mRay(r), mMinT(_min) {
+		SetMaxT(_max);
+	}
+	
+	RayInfo(VssRay *r,
+										const float _min,
+										const short _max
+										):mRay(r), mMaxT(_max) {
+		SetMinT(_min);
+	}
+	
+	friend bool operator<(const RayInfo &a, const RayInfo &b) {
+		return a.mRay < b.mRay;
+	}
+  
+	
+	float ExtrapOrigin(const int axis) const {
+		return mRay->GetOrigin(axis) + GetMinT()*mRay->GetDir(axis);
+	}
+	
+	float ExtrapTermination(const int axis) const {
+		return mRay->GetOrigin(axis) + GetMaxT()*mRay->GetDir(axis);
+	}
+	
+#if USE_FIXEDPOINT_T
+	float GetMinT () const { return mMinT/(float)(FIXEDPOINT_ONE); }
+	float GetMaxT () const { return mMaxT/(float)(FIXEDPOINT_ONE); }
+	
+	void SetMinT (const float t) {
+		mMinT = (short) (t*(float)(FIXEDPOINT_ONE));
+	}
+	
+	void SetMaxT (const float t) {
+		mMaxT = (short) (t*(float)(FIXEDPOINT_ONE));
+		mMaxT++;
+		//      if (mMaxT!=0xFFFF)
+		//	mMaxT++;
+	}
+#else
+	float GetMinT () const { return mMinT; }
+	float GetMaxT () const { return mMaxT; }
+	
+	void SetMinT (const float t) { mMinT = t; }
+	void SetMaxT (const float t) { mMaxT = t; }
+#endif 
+
+
+  int ComputeRayIntersection(const float axis,
+														 const float position,
+														 float &t
+														 ) const {
+		
+		// intersect the ray with the plane
+    float denom = mRay->GetDir(axis);
+    
+    if (fabs(denom) < 1e-20)
+      //if (denom == 0.0f)
+      return (mRay->GetOrigin(axis) > position) ? 1 : -1;
+    
+    t = (position - mRay->GetOrigin(axis))/denom;
+
+    if (t < GetMinT())
+      return (denom > 0) ? 1 : -1;
+
+    if (t > GetMaxT())
+      return (denom > 0) ? -1 : 1;
+
+		return 0;
+	}
+
+
+};
+
+
+	typedef vector<RayInfo> RayInfoContainer;
+	
+  enum { EInterior, ELeaf };
+
+  /////////////////////////////////
+  // The actual data goes here
+  
+  // link to the parent
+  VspKdTreeInterior *parent;
+
+  enum {SPLIT_X=0, SPLIT_Y, SPLIT_Z, SPLIT_DIRX, SPLIT_DIRY, SPLIT_DIRZ};
+  
+  // splitting axis
+  char axis;
+	
+  // depth
+  unsigned char depth;
+  
+  //  short depth;
+  //
+  /////////////////////////////////
+  
+  inline VspKdTreeNode(VspKdTreeInterior *p);
+
+  
+  virtual ~VspKdTreeNode() {};
+  virtual int Type() const  = 0;
+  
+
+  bool IsLeaf() const { return axis == -1; }
+  
+  virtual void Print(ostream &s) const = 0;
+
+  virtual int GetAccessTime() {
+    return 0x7FFFFFF;
+  }
+
+  
+	
+};
+
+// --------------------------------------------------------------
+// KD-tree node - interior node
+// --------------------------------------------------------------
+class VspKdTreeInterior :
+  public VspKdTreeNode
+{
+public:
+  // plane in local modelling coordinates
+  float position;
+
+  // pointers to children
+  VspKdTreeNode *back, *front;
+
+  // the bbox of the node
+  AxisAlignedBox3 bbox;
+
+  // the bbox of the node
+  AxisAlignedBox3 dirBBox;
+  
+  // data for caching
+  long accesses;
+  long lastAccessTime;
+  
+  VspKdTreeInterior(VspKdTreeInterior *p):VspKdTreeNode(p),
+					back(NULL),
+					front(NULL),
+					accesses(0),
+					lastAccessTime(-1)
+  { }
+
+  virtual int GetAccessTime() {
+    return lastAccessTime;
+  }
+
+  void SetupChildLinks(VspKdTreeNode *b, VspKdTreeNode *f) {
+    back = b;
+    front = f;
+    b->parent = f->parent = this;
+  }
+
+  void ReplaceChildLink(VspKdTreeNode *oldChild, VspKdTreeNode *newChild) {
+    if (back == oldChild)
+      back = newChild;
+    else
+      front = newChild;
+  }
+
+  virtual int Type() const  { return EInterior; }
+  
+  virtual ~VspKdTreeInterior() {
+    if (back)
+      delete back;
+    if (front)
+      delete front;
+  }
+  
+  virtual void Print(ostream &s) const {
+    if (axis == 0)
+      s<<"x ";
+    else
+      if (axis == 1)
+	s<<"y ";
+      else
+	s<<"z ";
+    s<<position<<" ";
+    back->Print(s);
+    front->Print(s);
+  }
+
+  
+	
+  int ComputeRayIntersection(const RayInfo &rayData,
+														 float &t
+														 ) {
+		return rayData.ComputeRayIntersection(axis, position, t);
+  }
+
+};
+
+
+// --------------------------------------------------------------
+// KD-tree node - leaf node
+// --------------------------------------------------------------
+class VspKdTreeLeaf :
+  public VspKdTreeNode
+{
+private:
+	int mPvsSize;
+public:
+  static int mailID;
+  int mailbox;
+  
+  RayInfoContainer rays;
+	bool mValidPvs;
+	
+  VspKdTreeLeaf(VspKdTreeInterior *p,
+							const int nRays
+							):VspKdTreeNode(p), rays(), mPvsSize(0), mValidPvs(false) {
+    rays.reserve(nRays);
+  }
+  
+  virtual ~VspKdTreeLeaf() { }
+
+  virtual int Type() const  { return ELeaf; }
+
+  virtual void Print(ostream &s) const {
+    s<<endl<<"L: r="<<rays.size()<<endl;
+  };
+  
+  void AddRay(const RayInfo &data) {
+		mValidPvs = false;
+    rays.push_back(data);
+    data.mRay->Ref();
+  }
+	
+	int GetPvsSize() const {
+		return mPvsSize;
+	}
+	void SetPvsSize(const int s) {
+		mPvsSize = s;
+	}
+
+	void
+	UpdatePvsSize();
+
+  void Mail() { mailbox = mailID; }
+  static void NewMail() { mailID++; }
+  bool Mailed() const { return mailbox == mailID; }
+  
+  bool Mailed(const int mail) {
+    return mailbox >= mailID + mail;
+  }
+
+	float GetAvgRayContribution() const {
+		return GetPvsSize()/((float)rays.size() + Limits::Small);
+	}
+
+	float GetSqrRayContribution() const {
+		return sqr(GetPvsSize()/((float)rays.size() + Limits::Small));
+	}
+
+};
+
+// Inline functions
+inline
+VspKdTreeNode::VspKdTreeNode(VspKdTreeInterior *p):
+  parent(p), axis(-1), depth(p ? p->depth + 1 : 0) {}
+
+
+
+// ---------------------------------------------------------------
+// Main LSDS search class
+// ---------------------------------------------------------------
+class VspKdTree
+{
+  struct TraversalData
+  {  
+    VspKdTreeNode *node;
+    AxisAlignedBox3 bbox;
+    int depth;
+    float priority;
+    
+    TraversalData() {}
+
+    TraversalData(VspKdTreeNode *n, const float p):
+      node(n), priority(p)
+    {}
+
+    TraversalData(VspKdTreeNode *n,
+		   const AxisAlignedBox3 &b,
+		   const int d):
+      node(n), bbox(b), depth(d) {}
+    
+		
+    // comparator for the 
+    struct less_priority : public
+    binary_function<const TraversalData, const TraversalData, bool> {
+			
+      bool operator()(const TraversalData a, const TraversalData b) {
+				return a.priority < b.priority;
+      }
+      
+    };
+
+    //    ~TraversalData() {}
+    //    TraversalData(const TraversalData &s):node(s.node), bbox(s.bbox), depth(s.depth) {}
+    
+    friend bool operator<(const TraversalData &a,
+													const TraversalData &b) {
+      //      return a.node->queries.size() < b.node->queries.size();
+      VspKdTreeLeaf *leafa = (VspKdTreeLeaf *) a.node;
+      VspKdTreeLeaf *leafb = (VspKdTreeLeaf *) b.node;
+#if 0
+			return
+				leafa->rays.size()*a.bbox.GetVolume()
+				<
+				leafb->rays.size()*b.bbox.GetVolume();
+#endif
+#if 1
+			return
+				leafa->GetPvsSize()*a.bbox.GetVolume()
+				<
+				leafb->GetPvsSize()*b.bbox.GetVolume();
+#endif
+#if 0
+			return
+				leafa->GetPvsSize()
+				<
+				leafb->GetPvsSize();
+#endif
+#if 0
+			return
+				leafa->GetPvsSize()/(leafa->rays.size()+1)
+				>
+				leafb->GetPvsSize()/(leafb->rays.size()+1);
+#endif
+#if 0
+			return
+				leafa->GetPvsSize()*leafa->rays.size()
+				<
+				leafb->GetPvsSize()*leafb->rays.size();
+#endif
+    }
+  };
+  
+  // simplified data for ray traversal only...
+
+  struct RayTraversalData {
+    
+    VspKdTreeNode::RayInfo rayData;
+    VspKdTreeNode *node;
+    
+    RayTraversalData() {}
+    RayTraversalData(VspKdTreeNode *n,
+										 const VspKdTreeNode::RayInfo &data):
+      rayData(data), node(n) {}
+  };
+	
+public:
+  /////////////////////////////
+  // The core pointer
+  VspKdTreeNode *root;
+  
+  /////////////////////////////
+  // Basic properties
+
+  // total number of nodes of the tree
+  int nodes;
+  // axis aligned bounding box of the scene
+  AxisAlignedBox3 bbox;
+
+  // axis aligned bounding box of directions
+  AxisAlignedBox3 dirBBox;
+  
+  /////////////////////////////
+  // Construction parameters
+
+  // epsilon used for the construction
+  float epsilon;
+
+  // ratio between traversal and intersection costs
+  float ct_div_ci;
+  // max depth of the tree
+  int termMaxDepth;
+  // minimal ratio of the volume of the cell and the query volume
+  float termMinSize;
+
+	// minimal pvs per node to still get subdivided
+  int termMinPvs;
+
+	// minimal ray number per node to still get subdivided
+  int termMinRays;
+	
+  // maximal cost ration to subdivide a node
+  float termMaxCostRatio;
+	
+	// maximal contribution per ray to subdivide the node
+	float termMaxRayContribution;
+
+	
+  // randomized construction
+  bool randomize;
+
+  // type of the splitting to use fo rthe tree construction
+  enum {ESplitRegular, ESplitHeuristic };
+  int splitType;
+	
+  // maximal size of the box on which the refdir splitting can be performed
+  // (relative to the scene bbox
+  float refDirBoxMaxSize;
+  
+  // maximum alovable memory in MB
+  float maxTotalMemory;
+
+  // maximum alovable memory for static kd tree in MB
+  float maxStaticMemory;
+
+  // this is used during the construction depending
+  // on the type of the tree and queries...
+  float maxMemory;
+
+
+  // minimal acess time for collapse
+  int accessTimeThreshold;
+
+	// minimal depth at which to perform collapse
+  int minCollapseDepth;
+
+  
+  // reusable array of split candidates
+  vector<SortableEntry> *splitCandidates;
+  /////////////////////////////
+
+  VspKdStatistics stat;
+	
+  
+  VspKdTree();
+  virtual ~VspKdTree();
+
+  virtual void
+  Construct(
+						VssRayContainer &rays,
+						AxisAlignedBox3 *forcedBoundingBox = NULL
+						);
+	
+  // incemental construction
+  virtual void UpdateRays(VssRayContainer &remove,
+													VssRayContainer &add
+													);
+
+	virtual void AddRays(
+											 VssRayContainer &add
+											 )
+	{
+		VssRayContainer remove;
+		UpdateRays(remove, add);
+	}
+
+  
+	
+  VspKdTreeNode *
+  Locate(const Vector3 &v);
+	
+  VspKdTreeNode *
+  SubdivideNode(VspKdTreeLeaf *leaf,
+								const AxisAlignedBox3 &box,
+								AxisAlignedBox3 &backBox,
+								AxisAlignedBox3 &frontBox
+								);
+	
+  VspKdTreeNode *
+  Subdivide(const TraversalData &tdata);
+	
+  int
+  SelectPlane(VspKdTreeLeaf *leaf,
+							const AxisAlignedBox3 &box,
+							float &position,
+							int &raysBack,
+							int &raysFront,
+							int &pvsBack,
+							int &pvsFront
+							);
+
+  void
+  SortSplitCandidates(
+											VspKdTreeLeaf *node,
+											const int axis
+											);
+	
+  
+  // return memory usage in MB
+  float GetMemUsage() const {
+    return 
+      (sizeof(VspKdTree) +
+       stat.Leaves()*sizeof(VspKdTreeLeaf) +
+       stat.Interior()*sizeof(VspKdTreeInterior) +
+       stat.rayRefs*sizeof(VspKdTreeNode::RayInfo))/(1024.0f*1024.0f);
+  }
+	
+  float GetRayMemUsage() const {
+    return 
+      stat.rays*(sizeof(VssRay))/(1024.0f*1024.0f);
+  }
+  
+	float
+	BestCostRatioHeuristic(
+												 VspKdTreeLeaf *node,
+												 int &axis,
+												 float &position,
+												 int &raysBack,
+												 int &raysFront,
+												 int &pvsBack,
+												 int &pvsFront
+												 );
+
+	float
+	BestCostRatioRegular(
+											 VspKdTreeLeaf *node,
+											 int &axis,
+											 float &position,
+											 int &raysBack,
+											 int &raysFront,
+											 int &pvsBack,
+											 int &pvsFront
+
+											 );
+	
+	float
+	EvalCostRatio(
+								VspKdTreeLeaf *node,
+								const int axis,
+								const float position,
+								int &raysBack,
+								int &raysFront,
+								int &pvsBack,
+								int &pvsFront
+								);
+
+  AxisAlignedBox3 GetBBox(const VspKdTreeNode *node) {
+    if (node->parent == NULL)
+      return bbox;
+
+    if (!node->IsLeaf())
+      return ((VspKdTreeInterior *)node)->bbox;
+
+    if (node->parent->axis >= 3)
+      return node->parent->bbox;
+      
+    AxisAlignedBox3 box(node->parent->bbox);
+    if (node->parent->front == node)
+      box.SetMin(node->parent->axis, node->parent->position);
+    else
+      box.SetMax(node->parent->axis, node->parent->position);
+    return box;
+  }
+
+  AxisAlignedBox3 GetDirBBox(const VspKdTreeNode *node) {
+
+    if (node->parent == NULL)
+      return dirBBox;
+    
+    if (!node->IsLeaf() )
+      return ((VspKdTreeInterior *)node)->dirBBox;
+
+    if (node->parent->axis < 3)
+      return node->parent->dirBBox;
+    
+    AxisAlignedBox3 dBBox(node->parent->dirBBox);
+
+    if (node->parent->front == node)
+      dBBox.SetMin(node->parent->axis - 3, node->parent->position);
+    else
+      dBBox.SetMax(node->parent->axis - 3, node->parent->position);
+    return dBBox;
+  }
+  
+  int
+  ReleaseMemory(const int time);
+
+  int
+  CollapseSubtree(VspKdTreeNode *node, const int time);
+
+  void
+  CountAccess(VspKdTreeInterior *node, const long time) {
+    node->accesses++;
+    node->lastAccessTime = time;
+  }
+
+  VspKdTreeNode *
+  SubdivideLeaf(
+								VspKdTreeLeaf *leaf,
+								const float SAThreshold
+								);
+
+  void
+  RemoveRay(VssRay *ray,
+						vector<VspKdTreeLeaf *> *affectedLeaves,
+						const bool removeAllScheduledRays
+						);
+
+  void
+  AddRay(VssRay *ray);
+	
+  void
+  TraverseInternalNode(
+											 RayTraversalData &data,
+											 stack<RayTraversalData> &tstack);
+
+	void
+  EvaluateLeafStats(const TraversalData &data);
+
+
+	int
+	GetRootPvsSize() const {
+		return GetPvsSize(root, bbox);
+	}
+	
+	int
+	GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const;
+
+	void
+	GetRayContributionStatistics(
+															 float &minRayContribution,
+															 float &maxRayContribution,
+															 float &avgRayContribution
+															 );
+
+	int
+	GenerateRays(const float ratioPerLeaf,
+							 SimpleRayContainer &rays);
+
+	float
+	GetAvgPvsSize();
+
+};
+
+
+#endif // __LSDS_KDTREE_H__
+
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.h	(revision 406)
+++ 	(revision )
@@ -1,513 +1,0 @@
-#ifndef _ViewCellKdTree_H__
-#define _ViewCellKdTree_H__
-
-#include <functional>
-using namespace std;
-
-#include "Containers.h"
-#include "AxisAlignedBox3.h"
-#include "Ray.h"
-#include "Pvs.h"
-
-  
-class ViewCellKdNode;
-class ViewCellKdLeaf;
-class ViewCellKdInterior;
-class Intersectable;
-
-
-// --------------------------------------------------------------
-// Static statistics for kd-tree search
-// --------------------------------------------------------------
-class ViewCellKdTreeStatistics 
-{
-public:  
-  // total number of nodes
-  int nodes;
-  // number of splits along each of the axes
-  int splits[7];
-  // totals number of rays
-  int rays;
-  // total number of query domains
-  int queryDomains;
-  // total number of ray references
-  int rayRefs;
-  // refs in non empty leafs
-  int rayRefsNonZeroQuery;
-  // total number of query references
-  int objectRefs;
-  // nodes with zero queries
-  int zeroQueryNodes;
-  // max depth nodes
-  int maxDepthNodes;
-  // max depth nodes
-  int minCostNodes;
-  // max number of rays per node
-  int maxObjectRefs;
-  // number of dynamically added ray refs
-  int addedRayRefs;
-  // number of dynamically removed ray refs
-  int removedRayRefs;
-  
-  // Constructor
-  ViewCellKdTreeStatistics() {
-    Reset();
-  }
-
-  int Nodes() const {return nodes;}
-  int Interior() const { return nodes/2; }
-  int Leaves() const { return (nodes/2) + 1; }
-
-  void Reset() {
-    nodes = 0;
-    for (int i=0; i<7; i++)
-      splits[i] = 0;
-    rays = queryDomains = 0;
-    rayRefs = rayRefsNonZeroQuery = objectRefs = 0;
-    zeroQueryNodes = 0;
-    maxDepthNodes = 0;
-    minCostNodes = 0;
-    maxObjectRefs = 0;
-    addedRayRefs = removedRayRefs = 0;
-  }
-
-  void
-  Print(ostream &app) const;
-
-  friend ostream &operator<<(ostream &s, const ViewCellKdTreeStatistics &stat) {
-    stat.Print(s);
-    return s;
-  }
-  
-};
-
-  
-class ViewCellKdInterior;
-/** Abstract class for kd-tree node */
-class ViewCellKdNode {
-public:
-  static int mailID;
-  int mailbox;
-  
-  void Mail() { mailbox = mailID; }
-  static void NewMail() { mailID++; }
-  bool Mailed() const { return mailbox == mailID; }
-
-
-  ViewCellKdNode(ViewCellKdInterior *parent);
-
-  /** Determines whether this node is a leaf or interior node
-      @return true if leaf
-  */
-  virtual bool IsLeaf() const = 0;
-
-  /** Determines whether this node is the root of the tree
-      @return true if root
-  */
-  virtual bool IsRoot() const {
-    return mParent == NULL;
-  }
-
-  /** Parent of the node - the parent is a little overhead for maintanance of the tree,
-      but allows various optimizations of tree traversal algorithms */
-  ViewCellKdInterior *mParent;
-  int mDepth;
-};
-
-/** Implementation of the kd-tree interior node */
-class ViewCellKdInterior : public ViewCellKdNode {
-
-public:
-    
-  ViewCellKdInterior(ViewCellKdInterior *parent):ViewCellKdNode(parent), mBack(NULL), mFront(NULL) {}
-
-  /** \sa ViewCellKdNode::IsLeaf() */
-  virtual bool IsLeaf() const { return false; }
-
-  /** splitting axis */
-  int mAxis;
-  /** splitting position, absolute position within the bounding box of this node */
-  float mPosition;
-  /** bounding box of interior node */
-  AxisAlignedBox3 mBox;
-
-  /** back node */
-  ViewCellKdNode *mBack;
-  /** front node */
-  ViewCellKdNode *mFront;
-
-  void SetupChildLinks(ViewCellKdNode *b, ViewCellKdNode *f) {
-    mBack = b;
-    mFront = f;
-    b->mParent = f->mParent = this;
-  }
-  
-  void ReplaceChildLink(ViewCellKdNode *oldChild, ViewCellKdNode *newChild) {
-    if (mBack == oldChild)
-      mBack = newChild;
-    else
-      mFront = newChild;
-  }
-  
-  
-};
-  
-  
-/** Implementation of the kd-tree leaf node */
-class ViewCellKdLeaf : public ViewCellKdNode {
-public:
-  ViewCellKdLeaf(ViewCellKdInterior *parent, const int objects):ViewCellKdNode(parent) {
-    mObjects.reserve(objects);
-  }
-  
-  void AddPassingRay(const Ray &ray,
-										 const int contributions) {
-    mPassingRays.AddRay(ray, contributions);
-		//		Debug << "adding passing ray" << endl;
-  }
-
-	
-	void AddPassingRay2(const Ray &ray,
-											const int objects,
-											const int viewcells
-											) {
-    mPassingRays.AddRay2(ray, objects, viewcells);
-		//		Debug << "adding passing ray" << endl;
-  }
-
-  /** \sa ViewCellKdNode::IsLeaf() */
-  virtual bool IsLeaf() const { return true; }
-  
-
-  
-
-  /** pointers to occluders contained in this node */
-  ObjectContainer mObjects;
-  
-  /** pointers to viewcells contained in this node */
-  //  ViewCellContainer mViewCells;
-
-  /** Ray set description of the rays passing through this node */
-  PassingRaySet mPassingRays;
-	
-  /** PVS consisting of visible ViewCellKd nodes */
-  ViewCellKdPvs mKdPvs;
-	
-	/** PVS consisting of visible objects */
-  ViewCellPvs mPvs;
-};
-
-  
-
-/** ViewCellKd for indexing scene entities - occluders/occludees/viewcells */
-class ViewCellKd {
-
-protected:
-  struct TraversalData
-  {  
-    ViewCellKdNode *mNode;
-    AxisAlignedBox3 mBox;
-    int mDepth;
-    float mPriority;
-    
-    TraversalData() {}
-
-    TraversalData(ViewCellKdNode *n, const float p):
-      mNode(n), mPriority(p)
-    {}
-    
-    TraversalData(ViewCellKdNode *n,
-		   const AxisAlignedBox3 &b,
-		   const int d):
-      mNode(n), mBox(b), mDepth(d) {}
-    
-
-    bool operator<(
-		   const TraversalData &b) const {
-      ViewCellKdLeaf *leafa = (ViewCellKdLeaf *) mNode;
-      ViewCellKdLeaf *leafb = (ViewCellKdLeaf *) b.mNode;
-      return 
-	leafa->mObjects.size()*mBox.SurfaceArea()
-	<
-	leafb->mObjects.size()*b.mBox.SurfaceArea();
-    }
-
-
-    // comparator for the 
-    struct less_priority : public
-    binary_function<const TraversalData, const TraversalData, bool> {
-      
-      bool operator()(const TraversalData a, const TraversalData b) {
-				return a.mPriority < b.mPriority;
-      }
-      
-    };
-
-  };
-
-  
-
-public:
-
-  enum {SPLIT_OBJECT_MEDIAN,
-	SPLIT_SPATIAL_MEDIAN,
-	SPLIT_SAH};
-  
-  ViewCellKd();
-  
-    
-  /** Insert view cell into the tree */
-  virtual void InsertViewCell(ViewCell *viewCell) {
-    //      mRoot->mViewcells.push_back(viewCell);
-  }
-
-  virtual bool Construct();
-
-  /** Check whether subdivision criteria are met for the given subtree.
-      If not subdivide the leafs of the subtree. The criteria are specified in
-      the environment as well as the subdivision method. By default surface area
-      heuristics is used.
-	
-      @param subtree root of the subtree
-
-      @return true if subdivision was performed, false if subdivision criteria
-      were already met
-  */
-  virtual ViewCellKdNode *Subdivide(const TraversalData &tdata);
-
-  /** Get the root of the tree */
-  ViewCellKdNode *GetRoot() const {
-    return mRoot;
-  }
-
-
-  AxisAlignedBox3 GetBox() const { return mBox; }
-
-  int
-  CastRay(
-	  Ray &ray
-	  );
-
-  const ViewCellKdTreeStatistics &GetStatistics() const {
-    return mStat;
-  }
-
-  void
-  CollectObjects(ViewCellKdNode *n, ObjectContainer &objects);
-	
-  void
-  CollectLeaves(vector<ViewCellKdLeaf *> &leaves);
-	
-  AxisAlignedBox3 GetBox(const ViewCellKdNode *node) const {
-    ViewCellKdInterior *parent = node->mParent;
-    if (parent == NULL)
-      return mBox;
-    
-    if (!node->IsLeaf())
-      return ((ViewCellKdInterior *)node)->mBox;
-    
-    AxisAlignedBox3 box(parent->mBox);
-    if (parent->mFront == node)
-      box.SetMin(parent->mAxis, parent->mPosition);
-    else
-      box.SetMax(parent->mAxis, parent->mPosition);
-    return box;
-  }
-	
-  ViewCellKdNode *
-  FindRandomNeighbor(ViewCellKdNode *n,
-										 bool onlyUnmailed
-										 );
-  
-  ViewCellKdNode *
-  ViewCellKd::GetRandomLeaf(const Plane3 &halfspace);
-
-	ViewCellKdNode *
-	GetRandomLeaf(const bool onlyUnmailed = false);
-
-  int
-  FindNeighbors(ViewCellKdNode *n,
-		vector<ViewCellKdNode *> &neighbors,
-		bool onlyUnmailed
-		);
-
-  int
-  CollectLeafPvs();
-
-protected:
-
-  struct RayData {
-    // pointer to the actual ray
-    Ray *ray;
-    
-    // endpoints  - do we need them?
-#if USE_FIXEDPOINT_T
-    short tmin, tmax;
-#else
-    float tmin, tmax;
-#endif
-
-    RayData():ray(NULL) {}
-    RayData(Ray *r):ray(r), tmin(0),
-
-#if USE_FIXEDPOINT_T
-#define FIXEDPOINT_ONE 0x7FFE
-			  //			  tmax(0xFFFF)
-			  tmax(FIXEDPOINT_ONE)
-#else
-      tmax(1.0f)
-#endif
-    {}
-
-    RayData(Ray *r,
-	    const float _min,
-	    const float _max
-	    ):ray(r) {
-      SetTMin(_min);
-      SetTMax(_max);
-    }
-
-    RayData(Ray *r,
-	    const short _min,
-	    const float _max
-	    ):ray(r), tmin(_min) {
-      SetTMax(_max);
-    }
-
-    RayData(Ray *r,
-	    const float _min,
-	    const short _max
-	    ):ray(r), tmax(_max) {
-      SetTMin(_min);
-    }
-
-    friend bool operator<(const RayData &a, const RayData &b) {
-      return a.ray < b.ray;
-    }
-    
-    
-    float ExtrapOrigin(const int axis) const {
-      return ray->GetLoc(axis) + GetTMin()*ray->GetDir(axis);
-    }
-    
-    float ExtrapTermination(const int axis) const {
-      return ray->GetLoc(axis) + GetTMax()*ray->GetDir(axis);
-    }
-    
-#if USE_FIXEDPOINT_T
-    float GetTMin () const { return tmin/(float)(FIXEDPOINT_ONE); }
-    float GetTMax () const { return tmax/(float)(FIXEDPOINT_ONE); }
-
-    void SetTMin (const float t) {
-      tmin = (short) (t*(float)(FIXEDPOINT_ONE));
-    }
-    
-    void SetTMax (const float t) {
-      tmax = (short) (t*(float)(FIXEDPOINT_ONE));
-      tmax++;
-      //      if (tmax!=0xFFFF)
-      //	tmax++;
-    }
-#else
-    float GetTMin () const { return tmin; }
-    float GetTMax () const { return tmax; }
-
-    void SetTMin (const float t) { tmin = t; }
-    void SetTMax (const float t) { tmax = t; }
-#endif 
-  };
-
-  struct RayTraversalData {
-    ViewCellKdNode *mNode;
-    Vector3 mExitPoint;
-    float mMaxT;
-    
-    RayTraversalData() {}
-    RayTraversalData(ViewCellKdNode *n,
-		     const Vector3 &p,
-		     const float maxt):
-      mNode(n), mExitPoint(p), mMaxT(maxt) {}
-  };
-
-  // --------------------------------------------------------------
-  // For sorting objects
-  // --------------------------------------------------------------
-  struct  SortableEntry
-  {
-    enum {
-      BOX_MIN,
-      BOX_MAX
-    };
-    
-    int type;
-    float value;
-    Intersectable *intersectable;
-    
-    SortableEntry() {}
-    SortableEntry(const int t, const float v, Intersectable *i):type(t),
-								value(v),
-								intersectable(i) {}
-    
-    bool operator<(const SortableEntry &b) const {
-      return value < b.value;
-    }
-    
-  };
-  
-  // reusable array of split candidates
-  vector<SortableEntry> *splitCandidates;
-
-  float
-  BestCostRatio(
-		ViewCellKdLeaf *node,
-		const AxisAlignedBox3 &box,
-		const int axis,
-		float &position,
-		int &objectsBack,
-		int &objectsFront
-		);
-
-  void
-  SortSplitCandidates(
-		      ViewCellKdLeaf *node,
-		      const int axis
-		      );
-
-  void
-  EvaluateLeafStats(const TraversalData &data);
-
-  ViewCellKdNode *
-  SubdivideNode(
-		ViewCellKdLeaf *leaf,
-		const AxisAlignedBox3 &box,
-		AxisAlignedBox3 &backBBox,
-		AxisAlignedBox3 &frontBBox
-		);
-
-  bool
-  TerminationCriteriaMet(const ViewCellKdLeaf *leaf);
-  
-  int
-  SelectPlane(ViewCellKdLeaf *leaf,
-	      const AxisAlignedBox3 &box,
-	      float &position
-	      );
-
-
-
-  float mSplitBorder;
-  int mTermMaxDepth;
-  int mTermMinCost;
-  float mMaxCostRatio;
-  float mCt_div_ci;
-  int mSplitMethod;
-  bool mSahUseFaces;
-  /// root of the tree
-  ViewCellKdNode *mRoot;
-  /// bounding box of the tree root
-  AxisAlignedBox3 mBox;
-  ViewCellKdTreeStatistics mStat;
-  
-};
-
-#endif
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.cpp	(revision 406)
+++ 	(revision )
@@ -1,840 +1,0 @@
-#include <stack>
-#include <algorithm>
-#include <queue>
-#include "Environment.h"
-#include "Mesh.h"
-#include "ViewCellKd.h"
-
-
-int ViewCellKdNode::mailID = 1;
-
-
-ViewCellKdNode::ViewCellKdNode(ViewCellKdInterior *parent):mParent(parent), mailbox(0)
-{
-  if (parent)
-    mDepth = parent->mDepth+1;
-  else
-    mDepth = 0;
-}
-
-ViewCellKd::ViewCellKd()
-{
-  mRoot = new ViewCellKdLeaf(NULL, 0);
-  environment->GetIntValue("ViewCellKd.Termination.maxDepth", mTermMaxDepth);
-  environment->GetIntValue("ViewCellKd.Termination.minCost", mTermMinCost);
-  environment->GetFloatValue("ViewCellKd.Termination.maxCostRatio", mMaxCostRatio);
-  environment->GetFloatValue("ViewCellKd.Termination.ct_div_ci", mCt_div_ci);
-  environment->GetFloatValue("ViewCellKd.splitBorder", mSplitBorder);
-
-  environment->GetBoolValue("ViewCellKd.sahUseFaces", mSahUseFaces);
-
-  char splitType[64];
-  environment->GetStringValue("ViewCellKd.splitMethod", splitType);
-  
-  mSplitMethod = SPLIT_SPATIAL_MEDIAN;
-  if (strcmp(splitType, "spatialMedian") == 0)
-    mSplitMethod = SPLIT_SPATIAL_MEDIAN;
-  else
-    if (strcmp(splitType, "objectMedian") == 0)
-      mSplitMethod = SPLIT_OBJECT_MEDIAN;
-    else
-      if (strcmp(splitType, "SAH") == 0)
-	mSplitMethod = SPLIT_SAH;
-      else {
-	cerr<<"Wrong kd split type "<<splitType<<endl;
-	exit(1);
-      }
-  splitCandidates = NULL;
-}
-
-bool
-ViewCellKd::Construct()
-{
-  if (!splitCandidates)
-    splitCandidates = new vector<SortableEntry>;
-
-  // first construct a leaf that will get subdivide
-  ViewCellKdLeaf *leaf = (ViewCellKdLeaf *) mRoot;
-
-  mStat.nodes = 1;
-
-  mBox.Initialize();
-  
-  ObjectContainer::const_iterator mi;
-  for ( mi = leaf->mObjects.begin();
-	mi != leaf->mObjects.end();
-	mi++) {
-    mBox.Include((*mi)->GetBox());
-  }
-
-  cout <<"ViewCellKd Root Box:"<< mBox<<endl;
-  mRoot = Subdivide(TraversalData(leaf, mBox, 0));
-
-  // remove the allocated array
-  delete splitCandidates;
-
-  return true;
-}
-
-ViewCellKdNode *
-ViewCellKd::Subdivide(const TraversalData &tdata)
-{
-
-  ViewCellKdNode *result = NULL;
-
-  priority_queue<TraversalData> tStack;
-  //  stack<STraversalData> tStack;
-  
-  tStack.push(tdata);
-  AxisAlignedBox3 backBox, frontBox;
-
-  
-  while (!tStack.empty()) {
-
-#if 0
-    if ( GetMemUsage() > maxMemory ) {
-      // count statistics on unprocessed leafs
-      while (!tStack.empty()) {
-	EvaluateLeafStats(tStack.top());
-	tStack.pop();
-      }
-      break;
-    }
-#endif
-
-    TraversalData data = tStack.top();
-    tStack.pop();
-    
-    ViewCellKdNode *node = SubdivideNode((ViewCellKdLeaf *) data.mNode,
-				 data.mBox,
-				 backBox,
-				 frontBox
-				 );
-    if (result == NULL)
-      result = node;
-    
-    if (!node->IsLeaf()) {
-
-      ViewCellKdInterior *interior = (ViewCellKdInterior *) node;
-      // push the children on the stack
-      tStack.push(TraversalData(interior->mBack, backBox, data.mDepth+1));
-      tStack.push(TraversalData(interior->mFront, frontBox, data.mDepth+1));
-      
-    } else {
-      EvaluateLeafStats(data);
-    }
-  }
-  
-  return result;
-
-}
-
-
-
-bool
-ViewCellKd::TerminationCriteriaMet(const ViewCellKdLeaf *leaf)
-{
-  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
-  return
-    (leaf->mObjects.size() <= mTermMinCost) ||
-    (leaf->mDepth >= mTermMaxDepth);
-  
-}
-
-
-int
-ViewCellKd::SelectPlane(ViewCellKdLeaf *leaf,
-		    const AxisAlignedBox3 &box,
-		    float &position
-		    )
-{
-  int axis = -1;
-  
-  switch (mSplitMethod)
-    {
-    case SPLIT_SPATIAL_MEDIAN: {
-      axis = box.Size().DrivingAxis();
-      position = (box.Min()[axis] + box.Max()[axis])*0.5f;
-      break;
-    }
-    case SPLIT_SAH: {
-      int objectsBack, objectsFront;
-      float costRatio;
-      bool mOnlyDrivingAxis = false;
-      if (mOnlyDrivingAxis) {
-				axis = box.Size().DrivingAxis();
-				costRatio = BestCostRatio(leaf,
-																	box,
-																	axis,
-																	position,
-																	objectsBack,
-																	objectsFront);
-      } else {
-				costRatio = MAX_FLOAT;
-				for (int i=0; i < 3; i++) {
-					float p;
-					float r = BestCostRatio(leaf,
-																	box,
-																	i,
-																	p,
-																	objectsBack,
-																	objectsFront);
-					if (r < costRatio) {
-						costRatio = r;
-						axis = i;
-						position = p;
-					}
-				}
-      }
-      
-      if (costRatio > mMaxCostRatio) {
-				//	cout<<"Too big cost ratio "<<costRatio<<endl;
-				axis = -1;
-      }
-      break;
-    }
-    
-    }
-  return axis;
-}
-
-ViewCellKdNode *
-ViewCellKd::SubdivideNode(
-		      ViewCellKdLeaf *leaf,
-		      const AxisAlignedBox3 &box,
-		      AxisAlignedBox3 &backBBox,
-		      AxisAlignedBox3 &frontBBox
-		      )
-{
-  
-  if (TerminationCriteriaMet(leaf))
-    return leaf;
-  
-  float position;
-  
-  // select subdivision axis
-  int axis = SelectPlane( leaf, box, position );
-
-  if (axis == -1) {
-    return leaf;
-  }
-  
-  mStat.nodes+=2;
-  mStat.splits[axis]++;
-  
-  // add the new nodes to the tree
-  ViewCellKdInterior *node = new ViewCellKdInterior(leaf->mParent);
-
-  node->mAxis = axis;
-  node->mPosition = position;
-  node->mBox = box;
-  
-  backBBox = box;
-  frontBBox = box;
-  
-  // first count ray sides
-  int objectsBack = 0;
-  int objectsFront = 0;
-  
-  backBBox.SetMax(axis, position);
-  frontBBox.SetMin(axis, position);
-
-  ObjectContainer::const_iterator mi;
-
-  for ( mi = leaf->mObjects.begin();
-	mi != leaf->mObjects.end();
-	mi++) {
-    // determine the side of this ray with respect to the plane
-    AxisAlignedBox3 box = (*mi)->GetBox();
-    if (box.Max(axis) > position ) 
-      objectsFront++;
-    
-    if (box.Min(axis) < position )
-      objectsBack++;
-  }
-
-  
-  ViewCellKdLeaf *back = new ViewCellKdLeaf(node, objectsBack);
-  ViewCellKdLeaf *front = new ViewCellKdLeaf(node, objectsFront);
-
-  // replace a link from node's parent
-  if (  leaf->mParent )
-    leaf->mParent->ReplaceChildLink(leaf, node);
-
-  // and setup child links
-  node->SetupChildLinks(back, front);
-  
-  for (mi = leaf->mObjects.begin();
-       mi != leaf->mObjects.end();
-       mi++) {
-    // determine the side of this ray with respect to the plane
-    AxisAlignedBox3 box = (*mi)->GetBox();
-
-    if (box.Max(axis) >= position )
-      front->mObjects.push_back(*mi);
-    
-    if (box.Min(axis) < position )
-      back->mObjects.push_back(*mi);
-    
-    mStat.objectRefs -= (int)leaf->mObjects.size();
-    mStat.objectRefs += objectsBack + objectsFront;
-  }
-  
-  delete leaf;
-  return node;
-}
-
-
-
-void
-ViewCellKdTreeStatistics::Print(ostream &app) const
-{
-  app << "===== ViewCellKd statistics ===============\n";
-
-  app << "#N_RAYS Number of rays )\n"
-      << rays <<endl;
-  app << "#N_DOMAINS  ( Number of query domains )\n"
-      << queryDomains <<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_MAXOBJECTREFS  ( Max number of rayRefs / leaf )\n" <<
-    maxObjectRefs << "\n";
-
-  app << "#N_NONEMPTYRAYREFS  ( Number of rayRefs in nonEmpty leaves / non empty leaf )\n" <<
-    rayRefsNonZeroQuery/(double)(Leaves() - zeroQueryNodes) << "\n";
-
-  app << "#N_LEAFDOMAINREFS  ( Number of query domain Refs / leaf )\n" <<
-    objectRefs/(double)Leaves() << "\n";
-
-  //  app << setprecision(4);
-
-  app << "#N_PEMPTYLEAVES  ( Percentage of leaves with zero query domains )\n"<<
-    zeroQueryNodes*100/(double)Leaves()<<endl;
-
-  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
-    maxDepthNodes*100/(double)Leaves()<<endl;
-
-  app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
-    minCostNodes*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 ViewCellKd statistics ==========\n";
-
-}
-
-
-
-void
-ViewCellKd::EvaluateLeafStats(const TraversalData &data)
-{
-
-  // the node became a leaf -> evaluate stats for leafs
-  ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)data.mNode;
-
-  if (data.mDepth > mTermMaxDepth)
-    mStat.maxDepthNodes++;
-  
-  if ( (int)(leaf->mObjects.size()) < mTermMinCost)
-    mStat.minCostNodes++;
-  
-  
-  if ( (int)(leaf->mObjects.size()) > mStat.maxObjectRefs)
-    mStat.maxObjectRefs = (int)leaf->mObjects.size();
-  
-}
-
-
-
-void
-ViewCellKd::SortSplitCandidates(
-			    ViewCellKdLeaf *node,
-			    const int axis
-			    )
-{
-  splitCandidates->clear();
-  
-  int requestedSize = 2 * (int)node->mObjects.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(ObjectContainer::const_iterator mi = node->mObjects.begin();
-      mi != node->mObjects.end();
-      mi++) {
-    AxisAlignedBox3 box = (*mi)->GetBox();
-
-    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MIN,
-					     box.Min(axis),
-					     *mi)
-			       );
-    
-    
-    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MAX,
-					     box.Max(axis),
-					     *mi)
-			       );
-  }
-  
-  stable_sort(splitCandidates->begin(), splitCandidates->end());
-}
-
-
-float
-ViewCellKd::BestCostRatio(
-											ViewCellKdLeaf *node,
-											const AxisAlignedBox3 &box,
-											const int axis,
-											float &position,
-											int &objectsBack,
-											int &objectsFront
-											)
-{
-
-  SortSplitCandidates(node, axis);
-  
-  // go through the lists, count the number of objects left and right
-  // and evaluate the following cost funcion:
-  // C = ct_div_ci  + (ol + or)/queries
-
-  float totalIntersections = 0.0f;
-  vector<SortableEntry>::const_iterator ci;
-
-  for(ci = splitCandidates->begin();
-      ci < splitCandidates->end();
-      ci++) 
-    if ((*ci).type == SortableEntry::BOX_MIN) {
-      totalIntersections += (*ci).intersectable->IntersectionComplexity();
-    }
-	
-  float intersectionsLeft = 0;
-  float intersectionsRight = totalIntersections;
-	
-  int objectsLeft = 0, objectsRight = (int)node->mObjects.size();
-  
-  float minBox = box.Min(axis);
-  float maxBox = box.Max(axis);
-  float boxArea = box.SurfaceArea();
-  
-  float minBand = minBox + mSplitBorder*(maxBox - minBox);
-  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
-  
-  float minSum = 1e20;
-  
-  for(ci = splitCandidates->begin();
-      ci < splitCandidates->end();
-      ci++) {
-    switch ((*ci).type) {
-    case SortableEntry::BOX_MIN:
-      objectsLeft++;
-      intersectionsLeft += (*ci).intersectable->IntersectionComplexity();
-      break;
-    case SortableEntry::BOX_MAX:
-      objectsRight--;
-      intersectionsRight -= (*ci).intersectable->IntersectionComplexity();
-      break;
-    }
-
-    if ((*ci).value > minBand && (*ci).value < maxBand) {
-      AxisAlignedBox3 lbox = box;
-      AxisAlignedBox3 rbox = box;
-      lbox.SetMax(axis, (*ci).value);
-      rbox.SetMin(axis, (*ci).value);
-
-      float sum;
-      if (mSahUseFaces)
-				sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
-      else
-				sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
-      
-      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
-      //      cout<<"cost= "<<sum<<endl;
-      
-      if (sum < minSum) {
-				minSum = sum;
-				position = (*ci).value;
-				
-				objectsBack = objectsLeft;
-				objectsFront = objectsRight;
-      }
-    }
-  }
-  
-  float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
-  float newCost = mCt_div_ci + minSum/boxArea;
-  float ratio = newCost/oldCost;
-  
-#if 0
-  cout<<"===================="<<endl;
-  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
-      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
-#endif
-  return ratio;
-}
-
-int
-ViewCellKd::CastRay(
-								Ray &ray
-								)
-{
-  int hits = 0;
-  
-  stack<RayTraversalData> tStack;
-  
-  float maxt = 1e6;
-  float mint = 0;
-
-  Intersectable::NewMail();
-
-  if (!mBox.GetMinMaxT(ray, &mint, &maxt))
-    return 0;
-
-  if (mint < 0)
-    mint = 0;
-  
-  maxt += Limits::Threshold;
-  
-  Vector3 entp = ray.Extrap(mint);
-  Vector3 extp = ray.Extrap(maxt);
-  
-  ViewCellKdNode *node = mRoot;
-  ViewCellKdNode *farChild;
-  float position;
-  int axis;
-  
-  while (1) {
-    if (!node->IsLeaf()) {
-      ViewCellKdInterior *in = (ViewCellKdInterior *) node;
-      position = in->mPosition;
-      axis = in->mAxis;
-
-      if (entp[axis] <= position) {
-				if (extp[axis] <= position) {
-					node = in->mBack;
-					// cases N1,N2,N3,P5,Z2,Z3
-					continue;
-				} else {
-					// case N4
-					node = in->mBack;
-					farChild = in->mFront;
-				}
-      }
-      else {
-				if (position <= extp[axis]) {
-					node = in->mFront;
-					// cases P1,P2,P3,N5,Z1
-					continue;
-				} else {
-					node = in->mFront;
-					farChild = in->mBack;
-					// case P4
-				}
-			}
-      // $$ modification 3.5.2004 - hints from Kamil Ghais
-      // case N4 or P4
-      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
-      tStack.push(RayTraversalData(farChild, extp, maxt));
-      extp = ray.GetLoc() + ray.GetDir()*tdist;
-      maxt = tdist;
-		} else {
-			// compute intersection with all objects in this leaf
-			ViewCellKdLeaf *leaf = (ViewCellKdLeaf *) node;
-//			if (ray.mFlags & Ray::STORE_KDLEAVES)
-//				ray.kdLeaves.push_back(leaf);
-			
-			ObjectContainer::const_iterator mi;
-			for ( mi = leaf->mObjects.begin();
-						mi != leaf->mObjects.end();
-						mi++) {
-				Intersectable *object = *mi;
-				if (!object->Mailed() ) {
-					object->Mail();
-					if (ray.mFlags & Ray::STORE_TESTED_OBJECTS)
-						ray.testedObjects.push_back(object);
-					hits += object->CastRay(ray);
-				}
-			}
-			
-			if (hits && ray.GetType() == Ray::LOCAL_RAY)
-				if (ray.intersections[0].mT <= maxt)
-					break;
-			
-			// get the next node from the stack
-			if (tStack.empty())
-				break;
-			
-			entp = extp;
-			mint = maxt;
-			if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
-				break;
-			
-			RayTraversalData &s  = tStack.top();
-			node = s.mNode;
-			extp = s.mExitPoint;
-			maxt = s.mMaxT;
-			tStack.pop();
-		}
-  }
-  return hits;
-}
-
-void
-ViewCellKd::CollectObjects(ViewCellKdNode *n, ObjectContainer &objects)
-{
-  stack<ViewCellKdNode *> nodeStack;
-
-  nodeStack.push(n);
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)node;
-      for (int j=0; j < leaf->mObjects.size(); j++) {
-				Intersectable *object = leaf->mObjects[j];
-				if (!object->Mailed()) {
-					object->Mail();
-					objects.push_back(object);
-				}
-      }
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      nodeStack.push(interior->mFront);
-      nodeStack.push(interior->mBack);
-    }
-  }
-}
-
-// Find random neighbor which was not mailed
-ViewCellKdNode *
-ViewCellKd::FindRandomNeighbor(ViewCellKdNode *n,
-													 bool onlyUnmailed
-													 )
-{
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-
-  AxisAlignedBox3 box = GetBox(n);
-  int mask = rand();
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
-	return node;
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      if (interior->mPosition > box.Max(interior->mAxis))
-	nodeStack.push(interior->mBack);
-      else
-	if (interior->mPosition < box.Min(interior->mAxis))
-	  nodeStack.push(interior->mFront);
-	else {
-	  // random decision
-	  if (mask&1)
-	    nodeStack.push(interior->mBack);
-	  else
-	    nodeStack.push(interior->mFront);
-	  mask = mask>>1;
-	}
-    }
-  }
-  
-  return NULL;
-}
-
-int
-ViewCellKd::FindNeighbors(ViewCellKdNode *n,
-		      vector<ViewCellKdNode *> &neighbors,
-		      bool onlyUnmailed
-		      )
-{
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-
-  AxisAlignedBox3 box = GetBox(n);
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      if ( node != n && (!onlyUnmailed || !node->Mailed()) ) 
-	neighbors.push_back(node);
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      if (interior->mPosition > box.Max(interior->mAxis))
-				nodeStack.push(interior->mBack);
-      else
-				if (interior->mPosition < box.Min(interior->mAxis))
-					nodeStack.push(interior->mFront);
-				else {
-					// random decision
-					nodeStack.push(interior->mBack);
-					nodeStack.push(interior->mFront);
-				}
-    }
-  }
-  
-  return (int)neighbors.size();
-}
-
-// Find random neighbor which was not mailed
-ViewCellKdNode *
-ViewCellKd::GetRandomLeaf(const Plane3 &plane)
-{
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-  
-  int mask = rand();
-  
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      return node;
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      ViewCellKdNode *next;
-	if (GetBox(interior->mBack).Side(plane) < 0)
-	  next = interior->mFront;
-	else
-	  if (GetBox(interior->mFront).Side(plane) < 0)
-	    next = interior->mBack;
-	  else {
-	    // random decision
-	    if (mask&1)
-	      next = interior->mBack;
-	    else
-	      next = interior->mFront;
-	    mask = mask>>1;
-	  }
-	nodeStack.push(next);
-    }
-  }
-  
-  
-  return NULL;
-}
-
-void
-ViewCellKd::CollectLeaves(vector<ViewCellKdLeaf *> &leaves)
-{
-  stack<ViewCellKdNode *> nodeStack;
-  nodeStack.push(mRoot);
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)node;
-      leaves.push_back(leaf);
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      nodeStack.push(interior->mBack);
-      nodeStack.push(interior->mFront);
-    }
-  }
-}
-
-
-int
-ViewCellKd::CollectLeafPvs()
-{
-  int totalPvsSize = 0;
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-  
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)node;
-      for (int j=0; j < leaf->mObjects.size(); j++) {
-	Intersectable *object = leaf->mObjects[j];
-	if (!object->Mailed()) {
-	  object->Mail();
-	  // add this node to pvs of all nodes it can see
-	 /* ViewCellKdPvsMap::iterator ni = object->mKdPvs.mEntries.begin();
-	  for (; ni != object->mKdPvs.mEntries.end(); ni++) {
-	    ViewCellKdNode *node = (*ni).first;
-	    // $$ JB TEMPORARY solution -> should add object PVS or explictly computed
-	    // kd tree PVS
-	    if (leaf->mKdPvs.AddSample(node))
-	      totalPvsSize++;
-	  }*/
-	}
-      }
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      nodeStack.push(interior->mFront);
-      nodeStack.push(interior->mBack);
-    }
-  }
-
-  return totalPvsSize;
-}
-
-
-ViewCellKdNode *
-ViewCellKd::GetRandomLeaf(const bool onlyUnmailed)
-{
-	stack<ViewCellKdNode *> nodeStack;
-  nodeStack.push(mRoot);
-
-	int mask = rand();
-	
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      if ( (!onlyUnmailed || !node->Mailed()) )
-				return node;
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-			// random decision
-			if (mask&1)
-				nodeStack.push(interior->mBack);
-			else
-				nodeStack.push(interior->mFront);
-			mask = mask>>1;
-		}
-	}
-  return NULL;
-}
