Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Makefile
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Makefile	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Makefile	(revision 382)
@@ -1,5 +1,5 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (1.07a) (Qt 3.3.2) on: Thu Nov 03 11:34:08 2005
+# Generated by qmake (1.07a) (Qt 3.3.2) on: Sat Nov 05 19:00:35 2005
 # Project:  preprocessor.pro
 # Template: app
@@ -415,4 +415,5 @@
 main.obj: main.cpp  \
 		SamplingPreprocessor.h \
+		VssPreprocessor.h \
 		ExactPreprocessor.h \
 		Parser.h \
@@ -441,4 +442,5 @@
 		common.h \
 		Ray.h \
+		VssRay.h \
 		
 
@@ -639,4 +641,5 @@
 		Environment.h \
 		VssRay.h \
+		Intersectable.h \
 		AxisAlignedBox3.h \
 		Statistics.h \
@@ -646,4 +649,5 @@
 		Matrix4x4.h \
 		Plane3.h \
+		Pvs.h \
 		
 
@@ -657,4 +661,6 @@
 		Polygon3.h \
 		ViewCell.h \
+		VssRay.h \
+		VssTree.h \
 		Containers.h \
 		AxisAlignedBox3.h \
@@ -672,4 +678,5 @@
 		Material.h \
 		Exporter.h \
+		Statistics.h \
 		
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.h	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.h	(revision 382)
@@ -43,18 +43,20 @@
   // totals number of rays
   int rays;
-  // total number of query domains
+	// initial size of the pvs
+  int initialPvsSize;
+	// 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 queryDomainRefs;
-  // nodes with zero queries
-  int zeroQueryNodes;
-  // max depth nodes
+
+	// max depth nodes
   int maxDepthNodes;
   // max depth nodes
   int minCostNodes;
+	// max ray contribution nodes
+  int maxRayContribNodes;
+	// max depth nodes
+  int minSizeNodes;
+	
   // max number of rays per node
   int maxRayRefs;
@@ -78,10 +80,12 @@
       splits[i] = 0;
     rays = queryDomains = 0;
-    rayRefs = rayRefsNonZeroQuery = queryDomainRefs = 0;
-    zeroQueryNodes = 0;
+    rayRefs = 0;
     maxDepthNodes = 0;
     minCostNodes = 0;
     maxRayRefs = 0;
     addedRayRefs = removedRayRefs = 0;
+		initialPvsSize = 0;
+		maxRayContribNodes = 0;
+		minSizeNodes = 0;
   }
 
@@ -101,5 +105,5 @@
 // For sorting rays
 // --------------------------------------------------------------
-struct  SSortableEntry
+struct  SortableEntry
 {
   enum EType {
@@ -112,10 +116,10 @@
   void *data;
   
-  SSortableEntry() {}
-  SSortableEntry(const int t, const float v, void *d):type(t),
-																											value(v),
-						      data(d) {}
-	
-  friend bool operator<(const SSortableEntry &a, const SSortableEntry &b) {
+  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;
   }
@@ -215,4 +219,30 @@
 	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;
+	}
+
+
 };
 
@@ -329,24 +359,9 @@
 
   
+	
   int ComputeRayIntersection(const RayInfo &rayData,
 														 float &t
 														 ) {
-		
-    // intersect the ray with the plane
-    float denom = rayData.mRay->GetDir(axis);
-    
-    if (fabs(denom) < 1e-20)
-      //if (denom == 0.0f)
-      return (rayData.mRay->GetOrigin(axis) > position) ? 1 : -1;
-    
-    t = (position - rayData.mRay->GetOrigin(axis))/denom;
-
-    if (t < rayData.GetMinT())
-      return (denom > 0) ? 1 : -1;
-
-    if (t > rayData.GetMaxT())
-      return (denom > 0) ? -1 : 1;
-    
-    return 0;
+		return rayData.ComputeRayIntersection(axis, position, t);
   }
 
@@ -365,8 +380,10 @@
   
   RayInfoContainer rays;
-
-  
-  VssTreeLeaf(VssTreeInterior *p, const int n):VssTreeNode(p), rays() {
-    rays.reserve(n);
+	int mPvsSize;
+  
+  VssTreeLeaf(VssTreeInterior *p,
+							const int nRays
+							):VssTreeNode(p), rays(), mPvsSize(0) {
+    rays.reserve(nRays);
   }
   
@@ -392,4 +409,7 @@
   }
 
+	float GetAvgRayContribution() const {
+		return mPvsSize/(float)rays.size();
+	}
 };
 
@@ -490,17 +510,22 @@
   // max depth of the tree
   int termMaxDepth;
-  // minimal cost of the node to still get subdivided
-  int termMinCost;
   // minimal ratio of the volume of the cell and the query volume
   float termMinSize;
 
-  // minimal cost imporvement to subdivide a node
-  float maxCostRatio;
-  
+	// minimal pvs per node to still get subdivided
+  int termMinPvs;
+
+  // 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, ESplitVolume, ESplitQueries };
+  enum {ESplitRegular, ESplitVolume };
   int splitType;
 	
@@ -528,5 +553,5 @@
   
   // reusable array of split candidates
-  vector<SSortableEntry> *splitCandidates;
+  vector<SortableEntry> *splitCandidates;
   /////////////////////////////
 
@@ -540,5 +565,5 @@
   Construct(
 						VssRayContainer &rays,
-						const bool onlyStaticPart
+						AxisAlignedBox3 *forcedBoundingBox = NULL
 						);
 	
@@ -565,7 +590,11 @@
   int
   SelectPlane(VssTreeLeaf *leaf,
-	      const AxisAlignedBox3 &box,
-	      float &position
-	      );
+							const AxisAlignedBox3 &box,
+							float &position,
+							int &raysBack,
+							int &raysFront,
+							int &pvsBack,
+							int &pvsFront
+							);
 
   void
@@ -590,4 +619,23 @@
   }
   
+	float
+	BestCostRatio(
+								VssTreeLeaf *node,
+								int &axis,
+								float &position,
+								int &raysBack,
+								int &raysFront
+								);
+	
+	float
+	EvalCostRatio(
+								VssTreeLeaf *node,
+								const int axis,
+								const float position,
+								int &raysBack,
+								int &raysFront,
+								int &pvsBack,
+								int &pvsFront
+								);
 
   AxisAlignedBox3 GetBBox(const VssTreeNode *node) {
@@ -664,4 +712,8 @@
   EvaluateLeafStats(const TraversalData &data);
 
+	void
+	EvalLeafPvs(VssTreeLeaf *leaf);
+
+	
 };
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 379)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 382)
@@ -1289,5 +1289,31 @@
 								 "-preprocessor=",
 								 "sampling");
-	
+
+
+	RegisterOption("VssTree.maxDepth", optInt, "kd_depth=", "12");
+  RegisterOption("VssTree.minPvs", optInt, "kd_minpvs=", "32");
+	RegisterOption("VssTree.maxCostRatio", optFloat, "maxcost=", "0.9");
+	RegisterOption("VssTree.maxRayContribution", optFloat, "maxraycontrib=", "0.5");
+
+	RegisterOption("VssTree.epsilon", optFloat, "kd_eps=", "1e-6");
+  RegisterOption("VssTree.ct_div_ci", optFloat, "kd_ctdivci=", "1.0");
+  RegisterOption("VssTree.randomize", optBool, "randomize", "false");
+  RegisterOption("VssTree.splitType", optString, "split=", "queries");
+  RegisterOption("VssTree.numberOfEndPointDomains", optInt, "endpoints=", "10000");
+
+  RegisterOption("VssTree.minSize", optFloat, "minsize=", "0.001");
+
+  RegisterOption("VssTree.maxTotalMemory", optFloat, "mem=", "60.0");
+  RegisterOption("VssTree.maxStaticMemory", optFloat, "statmem=", "8.0");
+
+  RegisterOption("VssTree.queryType", optString, "qtype=", "static");
+
+  RegisterOption("VssTree.queryPosWeight", optFloat, "qposweight=", "0.0");
+  RegisterOption("VssTree.useRefDirSplits", optBool, "refdir", "false");
+  RegisterOption("VssTree.refDirAngle", optFloat, "refangle=", "10");
+  RegisterOption("VssTree.refDirBoxMaxSize", optFloat, "refboxsize=", "0.1");
+  RegisterOption("VssTree.accessTimeThreshold", optInt, "accesstime=", "1000");
+  RegisterOption("VssTree.minCollapseDepth", optInt, "colldepth=", "4");
+
 }
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp	(revision 382)
@@ -22,4 +22,5 @@
 #include "Environment.h"
 #include "VssRay.h"
+#include "Intersectable.h"
 
 
@@ -32,4 +33,35 @@
 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
@@ -37,10 +69,11 @@
 {
   environment->GetIntValue("VssTree.maxDepth", termMaxDepth);
-
-	environment->GetIntValue("VssTree.minCost", termMinCost);
+	environment->GetIntValue("VssTree.minPvs", termMinPvs);
+	environment->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution);
+  environment->GetFloatValue("VssTree.maxCostRatio", termMaxCostRatio);
 
   environment->GetFloatValue("VssTree.minSize", termMinSize);
   termMinSize = sqr(termMinSize);
-
+	
   environment->GetFloatValue("VssTree.refDirBoxMaxSize", refDirBoxMaxSize);
   refDirBoxMaxSize = sqr(refDirBoxMaxSize);
@@ -48,9 +81,8 @@
   environment->GetFloatValue("VssTree.epsilon", epsilon);
   environment->GetFloatValue("VssTree.ct_div_ci", ct_div_ci);
-
+	
   environment->GetFloatValue("VssTree.maxTotalMemory", maxTotalMemory);
   environment->GetFloatValue("VssTree.maxStaticMemory", maxStaticMemory);
   
-  environment->GetFloatValue("VssTree.maxCostRatio", maxCostRatio);
   
 
@@ -77,13 +109,10 @@
     splitType = ESplitRegular;
   else
-    if (name.compare("volume") == 0)
+    if (name.compare("heuristics") == 0)
       splitType = ESplitVolume;
-    else
-      if (name.compare("queries") == 0)
-				splitType = ESplitQueries;
-      else {
-				cerr<<"Invalid VssTree split type "<<name<<endl;
-				exit(1);
-      }
+		else {
+			cerr<<"Invalid VssTree split type "<<name<<endl;
+			exit(1);
+		}
 	
   environment->GetBoolValue("VssTree.randomize", randomize);
@@ -91,5 +120,5 @@
   root = NULL;
   
-  splitCandidates = new vector<SSortableEntry>;
+  splitCandidates = new vector<SortableEntry>;
 }
 
@@ -109,8 +138,9 @@
   app << "===== VssTree statistics ===============\n";
 
-  app << "#N_RAYS Number of rays )\n"
+  app << "#N_RAYS ( Number of rays )\n"
       << rays <<endl;
-  app << "#N_DOMAINS  ( Number of query domains )\n"
-      << queryDomains <<endl;
+
+	app << "#N_INITPVS ( Initial PVS size )\n"
+      << initialPvsSize <<endl;
   
   app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
@@ -135,14 +165,6 @@
     maxRayRefs << "\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" <<
-    queryDomainRefs/(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"<<
@@ -151,4 +173,10 @@
   app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
     minCostNodes*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"<<
@@ -173,14 +201,10 @@
 VssTree::Construct(
 									 VssRayContainer &rays,
-									 const bool onlyStaticPart
+									 AxisAlignedBox3 *forcedBoundingBox
 									 )
 {
   stat.Start();
   
-  if (onlyStaticPart)
-    maxMemory = maxStaticMemory;
-  else
-    maxMemory = maxTotalMemory;
-
+	maxMemory = maxStaticMemory;
 
   if (root)
@@ -194,5 +218,4 @@
   
   bbox.Initialize();
-
   dirBBox.Initialize();
   
@@ -201,13 +224,27 @@
       ri++) {
     leaf->AddRay(VssTreeNode::RayInfo(*ri));
+
     bbox.Include((*ri)->GetOrigin());
     bbox.Include((*ri)->GetTermination());
     
-    dirBBox.Include( (*ri)->GetNormalizedDir() );
-  }
-  
-  
+		
+		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();
-  
+	EvalLeafPvs(leaf);
+  stat.initialPvsSize = leaf->mPvsSize;
   // Subdivide();
   root = Subdivide(TraversalData(leaf, bbox, 0));
@@ -216,11 +253,35 @@
     // force realease of this vector
     delete splitCandidates;
-    splitCandidates = new vector<SSortableEntry>;
+    splitCandidates = new vector<SortableEntry>;
   }
   
   stat.Stop();
 
-  cout<<"##################################"<<endl;
-  
+	stat.Print(cout);
+	cout<<"#Total memory="<<GetMemUsage()<<endl;
+
+}
+
+void
+VssTree::EvalLeafPvs(VssTreeLeaf *leaf)
+{
+	Intersectable::NewMail();
+	int pvsSize = 0;
+	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+			ri != leaf->rays.end();
+			ri++)
+		if ((*ri).mRay->IsActive()) {
+			Intersectable *object = (*ri).mRay->mOriginObject;
+			if (object && !object->Mailed()) {
+				pvsSize++;
+				object->Mail();
+			}
+			object = (*ri).mRay->mTerminationObject;
+			if (object && !object->Mailed()) {
+				pvsSize++;
+				object->Mail();
+			}
+		}
+	leaf->mPvsSize = pvsSize;
 }
 
@@ -283,29 +344,322 @@
 										 VssTreeLeaf *leaf,
 										 const AxisAlignedBox3 &box,
-										 float &position
+										 float &position,
+										 int &raysBack,
+										 int &raysFront,
+										 int &pvsBack,
+										 int &pvsFront
 										 )
 {
 
   int minDirDepth = 6;
-
-  Vector3 size = box.Size();
-  
-  float pDirSubdivision = 0.0f;
-  if (leaf->depth >= minDirDepth && RandomValue(0, 1) < pDirSubdivision) {
-    AxisAlignedBox3 dirbox = GetDirBBox(leaf);
-    Vector3 size = dirbox.Size();
-    int axis = size.DrivingAxis();
-    position = (dirbox.Min()[axis] + dirbox.Max()[axis])*0.5f;
-    return axis + 3;
-  }
-  
-  int axis = size.DrivingAxis();
-
-  
+  int axis;
+	float costRatio;
+	
   if (splitType == ESplitRegular) {
-		position = (box.Min()[axis] + box.Max()[axis])*0.5f;
-  }
-  
-  return axis;
+		int sRaysBack, sRaysFront;
+		int dRaysBack, dRaysFront;
+		int sPvsBack, sPvsFront;
+		int dPvsBack, dPvsFront;
+		
+		
+		int sAxis = box.Size().DrivingAxis();
+		float sPosition = (box.Min()[sAxis] + box.Max()[sAxis])*0.5f;
+		float sCostRatio = EvalCostRatio(leaf,
+																		 sAxis,
+																		 sPosition,
+																		 sRaysBack,
+																		 sRaysFront,
+																		 sPvsBack,
+																		 sPvsFront
+																		 );
+		//		cout<<"srays back="<<sRaysBack<<" rays front="<<sRaysFront<<" pvs back="<<sPvsBack<<
+		//			" pvs front="<<sPvsFront<<endl;
+
+		AxisAlignedBox3 dirBox = GetDirBBox(leaf);
+		int dAxis = dirBox.Size().DrivingAxis();
+		float dPosition = (dirBox.Min()[dAxis] + dirBox.Max()[dAxis])*0.5f;
+		float dCostRatio = EvalCostRatio(leaf,
+																		 dAxis+3,
+																		 dPosition,
+																		 dRaysBack,
+																		 dRaysFront,
+																		 dPvsBack,
+																		 dPvsFront
+																		 );
+		
+		//		cout<<"drays back="<<dRaysBack<<" rays front="<<dRaysFront<<" pvs back="<<dPvsBack<<
+		//			" pvs front="<<dPvsFront<<endl;
+
+		if (sCostRatio < dCostRatio) {
+			costRatio = sCostRatio;
+			axis = sAxis;
+			position = sPosition;
+			raysBack = sRaysBack;
+			raysFront = sRaysFront;
+			pvsBack = sPvsBack;
+			pvsFront = sPvsFront;
+			
+		} else {
+			costRatio = dCostRatio;
+			axis = dAxis+3;
+			position = dPosition;
+			raysBack = dRaysBack;
+			raysFront = dRaysFront;
+			pvsBack = dPvsBack;
+			pvsFront = dPvsFront;
+		}
+		
+	} else {
+    if (splitType == ESplitVolume)
+      costRatio = BestCostRatio(leaf,
+																axis,
+																position,
+																raysBack,
+																raysFront);
+		else {
+			cerr<<"VssTree: 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
+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;
+
+	float newCost;
+
+	Intersectable::NewMail(3);
+	
+	// eval pvs size
+	int pvsSize = leaf->mPvsSize;
+
+	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()) {
+				
+				// 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->mOriginObject, side, pvsBack, pvsFront);
+				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(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->mOriginObject, side, pvsBack, pvsFront);
+				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
+VssTree::BestCostRatio(
+											 VssTreeLeaf *node,
+											 int &axis,
+											 float &position,
+											 int &raysBack,
+											 int &raysFront
+											 )
+{
+	AxisAlignedBox3 box = GetBBox(node);
+	AxisAlignedBox3 dirBox = GetDirBBox(node);
+
+	axis = box.Size().DrivingAxis();
+	
+	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  + (ql*rl + qr*rr)/queries
+
+  int rl=0, rr=node->rays.size();
+	
+  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;
+  
+  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
+      ci < splitCandidates->end();
+      ci++) {
+    
+    switch ((*ci).type) {
+    case SortableEntry::ERayMin:
+      rl++;
+      break;
+    case SortableEntry::ERayMax:
+      rr--;
+      break;
+    }
+		
+    if ((*ci).value > minBand && (*ci).value < maxBand) {
+      
+      sum = rl*((*ci).value - minBox) + rr*(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;
+				
+      }
+    }
+  }
+  
+  float oldCost = node->rays.size();
+  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
+VssTree::SortSplitCandidates(
+														 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++) {
+    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());
 }
 
@@ -318,14 +672,22 @@
   VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
 
-  if (data.depth > termMaxDepth)
+  if (data.depth >= termMaxDepth)
     stat.maxDepthNodes++;
   
-  if ( (int)(leaf->rays.size()) < termMinCost)
-    stat.minCostNodes++;
-  
-
-  if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
+	//  if ( (int)(leaf->rays.size()) < termMinCost)
+	//    stat.minCostNodes++;
+	if ( leaf->mPvsSize < termMinPvs)
+		stat.minCostNodes++;
+  
+	if (leaf->GetAvgRayContribution() > termMaxRayContribution )
+		stat.maxRayContribNodes++;
+	
+	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
+		stat.minSizeNodes++;
+	}
+
+	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
     stat.maxRayRefs = leaf->rays.size();
-  
+
 }
 
@@ -341,13 +703,32 @@
 {
   
-  if ( ((int)leaf->rays.size() < termMinCost) ||
-       (leaf->depth >= termMaxDepth) ) {
-    return leaf;
-  }
-  
-  float position;
-  
+  if ( (leaf->mPvsSize < termMinPvs) ||
+			 (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);
+  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) {
@@ -364,5 +745,4 @@
   node->position = position;
   node->bbox = box;
-
   node->dirBBox = GetDirBBox(leaf);
   
@@ -370,64 +750,8 @@
   frontBBox = box;
 
-  // first count ray sides
-  int raysBack = 0;
-  int raysFront = 0;
-  
-  if (axis <= VssTreeNode::SPLIT_Z) {
-    backBBox.SetMax(axis, position);
-    frontBBox.SetMin(axis, position);
-    
-    // this is the main ray classification loop!
-    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 = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
-				
-				(*ri).mRay->mSide = side;
-				
-				if (side <= 0)
-					raysBack++;
-				
-				if (side >= 0)
-					raysFront++;
-				
-#if 0
-	cout<<"-------------------"<<endl;
-	cout<<"axis = "<<(int)node->axis<<"\t position="<<node->position<<endl;
-	cout<<"origin="<<(*ri).ray->GetOrigin()<<"\t  term="<<(*ri).ray->GetTermination()<<endl;
-	cout<<"side = "<<side<<"\t t="<<t<<endl;
-	cout<<"-------------------"<<endl;
-#endif
-      }
-  } 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->GetNormalizedDir(axis - 3) > position)
-					side = 1;
-				else
-					side = -1;
-				
-				
-				(*ri).mRay->mSide = side;
-				
-				if (side <= 0)
-					raysBack++;
-				
-				if (side >= 0)
-					raysFront++;
-      }
-  }
-  
   VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
+	back->mPvsSize = pvsBack;
   VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
+	front->mPvsSize = pvsFront;
 
   // replace a link from node's parent
@@ -436,8 +760,10 @@
   // and setup child links
   node->SetupChildLinks(back, front);
-
-  
+	
   if (axis <= VssTreeNode::SPLIT_Z) {
-    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		backBBox.SetMax(axis, position);
+    frontBBox.SetMin(axis, position);
+		
+		for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
 				ri != leaf->rays.end();
 				ri++) {
@@ -446,6 +772,9 @@
 				// first unref ray from the former leaf
 				(*ri).mRay->Unref();
-				
-				if ((*ri).mRay->mSide == 0) {
+
+				// 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(VssTreeNode::RayInfo((*ri).mRay,
@@ -465,5 +794,5 @@
 					}
 				} else
-					if ((*ri).mRay->mSide == 1) 
+					if (side == 1) 
 						front->AddRay(*ri);
 					else
@@ -475,7 +804,4 @@
     // rays front/back
     
-#if DEBUG_DIR_SPLIT
-    cout<<"dir split, depth="<<(int)leaf->depth<<" front= "<<raysFront<<" back="<<raysBack<<endl;
-#endif
     
     for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
@@ -485,6 +811,12 @@
 				// first unref ray from the former leaf
 				(*ri).mRay->Unref();
-				
-				if ((*ri).mRay->mSide == 1)
+
+				int side;
+				if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
+					side = 1;
+				else
+					side = -1;
+				
+				if (side == 1)
 					front->AddRay(*ri);
 				else
@@ -573,5 +905,6 @@
   // check if we should perform a dynamic subdivision of the leaf
   if (
-      leaf->rays.size() > (unsigned)termMinCost &&
+			//      leaf->rays.size() > (unsigned)termMinCost &&
+			(leaf->mPvsSize >= termMinPvs) &&
       SqrMagnitude(leafBBox.Size()) > sizeThreshold) {
     
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h	(revision 376)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h	(revision 382)
@@ -22,5 +22,5 @@
 	
 	// side of the ray - used for the ray classification
-  char mSide;
+	//  char mSide;
 	
 	// computed t
@@ -98,4 +98,6 @@
   Vector3 GetNormalizedDir() const { return (mTermination - mOrigin)*mInvSize; }
 
+	float GetDirParametrization(const int axis) const;
+	
 	float GetInvSize() const { return  mInvSize; }
   float GetOrigin(const int axis) const { return mOrigin[axis]; }
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp	(revision 376)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp	(revision 382)
@@ -8,4 +8,5 @@
 #include "ViewCell.h"
 #include "VssRay.h"
+#include "VssTree.h"
 
 VssPreprocessor::VssPreprocessor():
@@ -48,5 +49,5 @@
 	Intersectable *objectA, *objectB;
 	Vector3 pointA, pointB;
-
+	float bsize = Magnitude(box.Size());
 	if (mKdTree->CastRay(ray)) {
 		objectA = ray.intersections[0].mObject;
@@ -57,5 +58,10 @@
 		float tmin, tmax;
 		box.ComputeMinMaxT(ray, &tmin, &tmax);
+		if (tmax > bsize) {
+			cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
+			cerr<<"ray"<<ray<<endl;
+		}
 		pointA = ray.Extrap(tmax);
+		
 	}
 	
@@ -69,4 +75,9 @@
 		float tmin, tmax;
 		box.ComputeMinMaxT(ray, &tmin, &tmax);
+		if (tmax > bsize) {
+			cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
+			cerr<<"ray"<<ray<<endl;
+		}
+		
 		pointB = ray.Extrap(tmax);
 	}
@@ -86,12 +97,14 @@
 
 Vector3
-VssPreprocessor::GetViewpoint()
-{
-	AxisAlignedBox3 box = mKdTree->GetBox();
-
+VssPreprocessor::GetViewpoint(AxisAlignedBox3 *viewSpaceBox)
+{
+	AxisAlignedBox3 box;
+	
+	if (viewSpaceBox)
+		box =*viewSpaceBox;
+	else
+		box = mKdTree->GetBox();
+	
 	// shrink the box in the y direction
-	Vector3 diff(0, -box.Size().y*0.4f, 0);
-	box.Enlarge(diff);
-	
 	return box.GetRandomPoint();
 }
@@ -100,5 +113,5 @@
 VssPreprocessor::GetDirection(const Vector3 &viewpoint)
 {
-	int i = RandomValue(0, mObjects.size());
+	int i = RandomValue(0, mObjects.size()-1);
 	Intersectable *object = mObjects[i];
 	Vector3 point, normal;
@@ -114,6 +127,6 @@
   
   mSceneGraph->CollectObjects(&mObjects);
-	
-
+	cout<<"#NUM_OBJECTS (Total numner of objects)\n"<<mObjects.size()<<endl;
+	
   long startTime = GetTime();
   
@@ -121,5 +134,13 @@
   int totalSamples = 0;
 
-
+	AxisAlignedBox3 *viewSpaceBox = NULL;
+
+	AxisAlignedBox3 box = mKdTree->GetBox();
+	box.Enlarge(box.Size()*-Vector3(0.45, 0.45, 0.45));
+
+	bool useViewSpaceBox = false;
+	if (useViewSpaceBox)
+		viewSpaceBox = &box;
+	
   while (totalSamples < mTotalSamples) {
 		int passContributingSamples = 0;
@@ -132,8 +153,8 @@
 		
 		for (int k=0; k < mSamplesPerPass; k++) {
-
-			Vector3 viewpoint = GetViewpoint();
+			
+			Vector3 viewpoint = GetViewpoint(viewSpaceBox);
 			Vector3 direction = GetDirection(viewpoint);
-
+			
 			VssRay *vssRay = CastRay(viewpoint, direction);
 			
@@ -153,5 +174,5 @@
     
 		mPass++;
-
+		
 		int pvsSize = 0;
 		float avgRayContrib = (passContributingSamples > 0) ? 
@@ -176,6 +197,10 @@
 	
 	cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
-	cout << "#totalRayStackSize=" << mVssRays.size() << endl;
-
+	cout << "#totalRayStackSize=" << mVssRays.size() << endl <<flush;
+
+	VssTree *vssTree = new VssTree;
+
+	vssTree->Construct(mVssRays, viewSpaceBox);
+	
   return true;
 }
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.cpp	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.cpp	(revision 382)
@@ -141,2 +141,16 @@
   return SqrDistance(point, center) < sqrRadius;
 }
+
+
+float
+VssRay::GetDirParametrization(const int axis) const
+{
+	float angle;
+	Vector3 dir = GetNormalizedDir();
+	if (axis == 0) {
+		angle = atan2(dir.x, dir.z);
+	} else {
+		angle = asin(dir.y);
+	}
+	return angle;
+}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h	(revision 354)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h	(revision 382)
@@ -11,4 +11,5 @@
 	// both mailId and mailbox should be unique for each thread!!!
 	static int sMailId;
+	static int sReservedMailboxes;
 	int mMailbox;
 	
@@ -26,9 +27,20 @@
 	int GetId() { return mId; }
 	
+
+	static void NewMail(const int reserve = 1) {
+		sMailId += sReservedMailboxes;
+		sReservedMailboxes = reserve;
+	}
+	
 	void Mail() { mMailbox = sMailId; }
-  static void NewMail() { sMailId++; }
   bool Mailed() const { return mMailbox == sMailId; }
-  int IncMail() { return ++mMailbox - sMailId; }
 
+	void Mail(const int mailbox) { mMailbox = sMailId + mailbox; }
+	bool Mailed(const int mailbox) const { return mMailbox == sMailId + mailbox; }
+
+	int IncMail() { return ++mMailbox - sMailId; }
+
+
+	
   virtual AxisAlignedBox3 GetBox() = 0;
   
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/default.env
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/default.env	(revision 376)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/default.env	(revision 382)
@@ -9,9 +9,9 @@
 #	filename vienna.x3d
 # filename ../data/vienna/vienna-simple.x3d
-filename ../data/vienna/vienna-buildings.x3d
+# filename ../data/vienna/vienna-buildings.x3d
 # filename ../data/vienna/viewcells-25-sel.x3d
 # filename ../data/atlanta/atlanta2.x3d
 #	filename ../data/soda/soda.dat
-# filename ../data/soda/soda5.dat
+filename ../data/soda/soda5.dat
 }
 
@@ -22,6 +22,28 @@
 
 VssPreprocessor {
-	totalSamples 1000000
-	samplesPerPass  100000
+	totalSamples 200000
+	samplesPerPass  50000
+}
+
+VssTree {
+	epsilon		1e-6
+
+	maxDepth 	40
+	minPvs		10
+	minSize		0.00001
+	maxCostRatio	1.0
+	maxRayContribution 0.5
+	
+	maxTotalMemory	500
+	maxStaticMemory	200
+
+	splitType regular
+#	splitType heuristics
+
+	numberOfEndPointDomains	10000
+	ct_div_ci	0.0
+	randomize	false
+
+	refDirBoxMaxSize	0.1
 }
 
@@ -66,5 +88,5 @@
 
 Sampling {
-	totalSamples 100000000
+	totalSamples 10000000
 	samplesPerPass  3
 }
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro	(revision 382)
@@ -40,3 +40,2 @@
 ViewCell.cpp ViewCellBsp.cpp Halton.cpp VssRay.cpp VssTree.cpp VssPreprocessor.cpp
 
-
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.h	(revision 376)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.h	(revision 382)
@@ -28,5 +28,5 @@
 
 	Vector3
-	GetViewpoint();
+	GetViewpoint(AxisAlignedBox3 *viewSpaceBox);
 
 	Vector3
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp	(revision 376)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp	(revision 382)
@@ -5,4 +5,5 @@
 
 int Intersectable::sMailId = 21843194198;
+int Intersectable::sReservedMailboxes = 1;
 
 void
