Index: trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj	(revision 390)
+++ trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj	(revision 406)
@@ -243,4 +243,10 @@
 			</File>
 			<File
+				RelativePath="..\src\RenderSimulator.cpp">
+			</File>
+			<File
+				RelativePath="..\src\RenderSimulator.h">
+			</File>
+			<File
 				RelativePath="..\src\SamplingPreprocessor.cpp">
 			</File>
@@ -283,4 +289,10 @@
 			<File
 				RelativePath="..\src\ViewCellBsp.h">
+			</File>
+			<File
+				RelativePath="..\src\ViewCellKd.cpp">
+			</File>
+			<File
+				RelativePath="..\src\ViewCellKd.h">
 			</File>
 			<File
@@ -358,7 +370,4 @@
 				RelativePath="..\include\ViewCell.h">
 			</File>
-			<File
-				RelativePath="..\include\ViewCellBsp.h">
-			</File>
 		</Filter>
 	</Files>
Index: trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 404)
+++ trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 406)
@@ -56,5 +56,5 @@
 
 Sampling {
-	totalSamples	300000
+	totalSamples	1000000
 	samplesPerPass  3
 }
@@ -126,12 +126,12 @@
 	#splitPlaneStrategy 130
 	
-	splitPlaneStrategy 1024
+	splitPlaneStrategy 1
 	
-	maxPolyCandidates 0
-	maxRayCandidates 50
+	maxPolyCandidates 50
+	maxRayCandidates 100
 	
 	Termination {
 		# autopartition
-		maxRays 300
+		maxRays 200
 		maxPolygons -1
 		maxDepth 40
@@ -161,2 +161,7 @@
 	}
 }
+
+Simulation {
+	objRenderCost 1.0
+	vcOverhead 0.05
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h	(revision 385)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h	(revision 406)
@@ -8,5 +8,5 @@
 class Ray;
 class Intersectable;
-
+class ViewCellKdNode;
 
 template<typename T>
@@ -142,8 +142,9 @@
 
 typedef std::map<KdNode *, PvsData<KdNode *>, LtSample<KdNode *> > KdPvsMap;
+typedef std::map<ViewCellKdNode *, PvsData<ViewCellKdNode *>, LtSample<ViewCellKdNode *> > ViewCellKdPvsMap;
 typedef std::map<Intersectable *, PvsData<Intersectable *>, LtSample<Intersectable *> > ViewCellPvsMap;
 typedef PvsData<Intersectable *> ViewCellPvsData;
 typedef PvsData<KdNode *> KdPvsData;
-//typedef Pvs<KdNode *> KdPvs;
+typedef Pvs<ViewCellKdNode *> ViewCellKdPvs;
 typedef Pvs<Intersectable *> ViewCellPvs;
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h	(revision 406)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h	(revision 406)
@@ -0,0 +1,56 @@
+#ifndef _RenderSimulator_h__
+#define _RenderSimulator_h__
+
+#include "common.h"
+
+class BspTree;
+class KdTree;
+class ViewCell;
+
+
+/** Simulated rendering using a simple render heuristics. Used to evaluate the 
+	quality of the view cell partition.
+*/
+class RenderSimulator
+{
+public:
+	RenderSimulator();
+	RenderSimulator(float objRendercost, float vcOverhead);
+	virtual Real SimulateRendering() = 0;
+
+protected:
+	// render time for single object of the PVS
+	float mObjRenderCost;
+	// const overhead for crossing a view cell border
+	float mVcOverhead;
+};
+
+class BspTreeRenderSimulator: public RenderSimulator
+{
+public:
+	BspTreeRenderSimulator(float objRenderCost, float vcOverhead, BspTree *bspTree);
+
+	Real SimulateRendering();
+
+protected:
+	/** Simulates rendering of the pvs of one view cell, with given rendering time for an object.
+		@param viewCell the view cell holding the Pvs
+		@param objRenderTime estimated render time for one object of the Pvs
+	*/
+	Real RenderPvs(ViewCell &viewCell, const float objRenderTime) const;
+
+	BspTree *mBspTree;
+};
+
+class KdTreeRenderSimulator: public RenderSimulator
+{
+public:
+	KdTreeRenderSimulator(float objRenderCost, float vcOverhead, KdTree *kdTree);
+
+	Real SimulateRendering();
+
+protected:
+	KdTree *mKdTree;
+};
+
+#endif // RenderSimulator
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 400)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 406)
@@ -814,4 +814,5 @@
 	static int sMinPvs;
 	static int sMaxPvs;
+	static bool sPvsUseArea;
 
 private:
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 401)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 406)
@@ -214,5 +214,9 @@
   /** Classifies ray with respect to the plane.
   */
-  int ClassifyPlane(const Plane3 &plane, const float minT, const float maxT) const;
+  int ClassifyPlane(const Plane3 &plane, 
+					const float minT, 
+					const float maxT,
+					Vector3 &entP,
+					Vector3 &extP) const;
 
 private:
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h	(revision 390)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h	(revision 406)
@@ -192,5 +192,5 @@
   virtual float IntersectionComplexity() {  return (float)mMesh->mFaces.size(); }
 
-	virtual int NumberOfFaces() const {  return mMesh->mFaces.size(); }
+	virtual int NumberOfFaces() const {  return (int)mMesh->mFaces.size(); }
 
   virtual int Type() const { return MESH_INSTANCE; }
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 403)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 406)
@@ -1177,15 +1177,15 @@
   RegisterOption("ViewCells.PostProcessing.minPvsDif",
 		 optInt,
-		 "-view_cells_min_pvs_dif",
+		 "-view_cells_post_processing_min_pvs_dif",
 		 "10");
-  
+
   RegisterOption("ViewCells.PostProcessing.maxPvs",
 		 optInt,
-		 "-view_cells_max_pvs",
+		 "-view_cells_post_processing_max_pvs",
 		 "300");
 
   RegisterOption("ViewCells.PostProcessing.minPvs",
 		 optString,
-		 "-view_cells_min_pvs",
+		 "-view_cells_post_processing_min_pvs",
 		 "10");
 
@@ -1199,4 +1199,14 @@
 		 "-view_cells_height=",
 		 "5.0");
+   
+  RegisterOption("Simulation.objRenderCost",
+		 optFloat,
+		 "-simulation_obj_render_cost",
+		 "1.0");
+
+  RegisterOption("Simulation.vcOverhead",
+		 optFloat,
+		 "-simulation_vc_overhead",
+		 "0.05");
 
   RegisterOption("BspTree.Construction.input",
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 404)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 406)
@@ -67,4 +67,6 @@
 */
 float BspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
+
+bool BspTree::sPvsUseArea = true;
 
 /****************************************************************/
@@ -1036,6 +1038,5 @@
 }
 
-Plane3 BspTree::SelectPlane(BspLeaf *leaf, 
-							BspTraversalData &data)
+Plane3 BspTree::SelectPlane(BspLeaf *leaf, BspTraversalData &data)
 {
 	if (data.mPolygons->empty() && data.mRays->empty())
@@ -1071,10 +1072,26 @@
 	if (sSplitPlaneStrategy & RANDOM_POLYGON)
 	{
-		if (!data.mPolygons->empty())
+        if (!data.mPolygons->empty())
 		{
 			Polygon3 *nextPoly = (*data.mPolygons)[Random((int)data.mPolygons->size())];
 			return nextPoly->GetSupportingPlane();
 		}
-		
+		else
+		{
+			const int candidateIdx = Random((int)data.mRays->size());
+			BoundedRay *bRay = (*data.mRays)[candidateIdx];
+
+			Ray *ray = bRay->mRay;
+						
+			const Vector3 minPt = ray->Extrap(bRay->mMinT);
+			const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
+
+			const Vector3 pt = (maxPt + minPt) * 0.5;
+
+			const Vector3 normal = ray->GetDir();
+			
+			return Plane3(normal, pt);
+		}
+
 		return Plane3();
 	}
@@ -1084,10 +1101,10 @@
 }
 
-Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf,
-									  BspTraversalData &data)
+Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf, BspTraversalData &data)
 {
 	float lowestCost = MAX_FLOAT;
 	Plane3 bestPlane;
-	
+	Plane3 plane;
+
 	int limit = Min((int)data.mPolygons->size(), sMaxPolyCandidates);
 	
@@ -1096,5 +1113,4 @@
 	for (int i = 0; i < limit; ++ i)
 	{
-		Debug << i << " ";
 		candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
 		
@@ -1116,57 +1132,24 @@
 	const BoundedRayContainer *rays = data.mRays;
 
-	for (int i = 0; i < sMaxRayCandidates; ++ i)
-	{
-		Plane3 plane;
-
-		if (0)
-		{
-			Vector3 pt[3];
-			int idx[3];
-			int cmaxT = 0;
-			int cminT = 0;
-			bool chooseMin = false;
-
-			for (int j = 0; j < 3; j ++)
-			{
-				idx[j] = Random((int)rays->size() * 2);
-				
-				if (idx[j] >= (int)rays->size())
-				{
-					idx[j] -= (int)rays->size();
-					
-					chooseMin = (cminT < 2);
-				}
-				else
-					chooseMin = (cmaxT < 2);
-
-				BoundedRay *bRay = (*rays)[idx[j]];
-				pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : bRay->mRay->Extrap(bRay->mMaxT);
-			}	
+	for (int i = 0; i < sMaxRayCandidates / 2; ++ i)
+	{
+		candidateIdx = Random((int)rays->size());
+		BoundedRay *bRay = (*rays)[candidateIdx];
+
+		Ray *ray = bRay->mRay;
+						
+		const Vector3 minPt = ray->Extrap(bRay->mMinT);
+		const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
+
+		const Vector3 pt = (maxPt + minPt) * 0.5;
+
+		const Vector3 normal = ray->GetDir();
 			
-			plane = Plane3(pt[0], pt[1], pt[2]);
-		}
-		else
-		{
-			candidateIdx = Random((int)rays->size());
-			BoundedRay *bRay = (*rays)[candidateIdx];
-
-			Ray *ray = bRay->mRay;
-						
-			const Vector3 minPt = ray->Extrap(bRay->mMinT);
-			const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
-
-			const Vector3 pt = (maxPt + minPt) * 0.5;
-
-			const Vector3 normal = ray->GetDir();
-			
-			plane = Plane3(normal, pt);
-		}
-		
+		plane = Plane3(normal, pt);
+	
 		const float candidateCost = SplitPlaneCost(plane, data);
 
 		if (candidateCost < lowestCost)
 		{
-			Debug << "choose ray plane: " << candidateCost << endl;
 			bestPlane = plane;
 			
@@ -1175,5 +1158,45 @@
 	}
 
-	Debug << "Plane lowest cost: " << lowestCost << endl;
+	Debug << "lowest: " << lowestCost << endl;
+		
+	for (int i = 0; i < sMaxRayCandidates / 2; ++ i)
+	{
+		Vector3 pt[3];
+		int idx[3];
+		int cmaxT = 0;
+		int cminT = 0;
+		bool chooseMin = false;
+
+		for (int j = 0; j < 3; j ++)
+		{
+			idx[j] = Random((int)rays->size() * 2);
+				
+			if (idx[j] >= (int)rays->size())
+			{
+				idx[j] -= (int)rays->size();
+				
+				chooseMin = (cminT < 2);
+			}
+			else
+				chooseMin = (cmaxT < 2);
+
+			BoundedRay *bRay = (*rays)[idx[j]];
+			pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : bRay->mRay->Extrap(bRay->mMaxT);
+		}	
+			
+		plane = Plane3(pt[0], pt[1], pt[2]);
+
+		const float candidateCost = SplitPlaneCost(plane, data);
+
+		if (candidateCost < lowestCost)
+		{
+			Debug << "choose ray plane 2: " << candidateCost << endl;
+			bestPlane = plane;
+			
+			lowestCost = candidateCost;
+		}
+	}	
+
+	Debug << "plane lowest cost: " << lowestCost << endl;
 	return bestPlane;
 }
@@ -1358,5 +1381,5 @@
 		Intersectable::NewMail(); frontAndBackId = ViewCell::sMailId;
 
-		if (1) // use front and back cell areas to approximate volume
+		if (sPvsUseArea) // use front and back cell areas to approximate volume
 		{	
 			// construct child geometry with regard to the candidate split plane
@@ -1371,8 +1394,4 @@
 			pOverall = area;
 		}
-		else
-		{
-			pOverall = (float)rays.size();
-		}
 	}
 			
@@ -1385,6 +1404,8 @@
 		const float maxT = (*rit)->mMaxT;
 
+		Vector3 entP, extP;
+
 		const int cf = 
-			ray->ClassifyPlane(candidatePlane, minT, maxT);
+			ray->ClassifyPlane(candidatePlane, minT, maxT, entP, extP);
 
 		if (sSplitPlaneStrategy & LEAST_RAY_SPLITS)
@@ -1416,11 +1437,52 @@
 			}
 
-			// use #rays to approximate volume
-			if (0)
+			if (!sPvsUseArea) // use front and back cell areas to approximate volume
 			{
-				if ((cf == Ray::BACK) || (cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT))
-					++ pFront;
-				if ((cf == Ray::FRONT) || (cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT))
-					++ pBack;
+				float len = Distance(entP, extP);
+				pOverall += len;
+
+				// use length of rays to approximate volume
+				switch (cf)
+				{
+					case Ray::COINCIDENT:
+						pBack += len;
+						pFront += len;						
+						break;
+					case Ray::BACK:
+						pBack += len;
+						break;
+					case Ray::FRONT:
+						pFront += len;
+						break;
+					case Ray::FRONT_BACK:
+						{
+							// find intersection of ray segment with plane
+							const Vector3 extp = ray->Extrap(maxT);
+							const float t = candidatePlane.FindT(ray->GetLoc(), extp);
+					
+							const float newT = t * maxT;
+							float newLen = Distance(ray->Extrap(newT), extp);
+
+							pFront += len - newLen;
+							pBack += newLen;
+						}
+						break;
+					case Ray::BACK_FRONT:
+						{
+							// find intersection of ray segment with plane
+							const Vector3 extp = ray->Extrap(maxT);
+							const float t = candidatePlane.FindT(ray->GetLoc(), extp);
+					
+							const float newT = t * maxT;
+							float newLen = Distance(ray->Extrap(newT), extp);
+
+							pFront += len;
+							pBack += len - newLen;
+						}
+						break;
+					default:
+						Debug << "Should not come here" << endl;
+						break;
+				}
 			}
 		}
@@ -1443,49 +1505,10 @@
 	}
 
-if ((val < 0.00001) || (frontPvs + backPvs < pvs))
-{
-	// create three unique ids for pvs heuristics
-	Intersectable::NewMail(); backId = ViewCell::sMailId;
-	Intersectable::NewMail(); frontId = ViewCell::sMailId;
-	Intersectable::NewMail(); frontAndBackId = ViewCell::sMailId;
-
-
-	BoundedRayContainer::const_iterator rit, rit_end = rays.end();
-
-	int fpvs = 0; int bpvs = 0;
-	for (rit = rays.begin(); rit != rays.end(); ++ rit)
-	{
-		Ray *ray = (*rit)->mRay;
-		const float minT = (*rit)->mMinT;
-		const float maxT = (*rit)->mMaxT;
-
-		const int cf = 
-			ray->ClassifyPlane(candidatePlane, minT, maxT);
-
-		Debug << "classification: " << cf << endl;
-		if (!ray->intersections.empty())
-		{
-			// in case the ray intersects an objcrs
-			// assure that we only count a object 
-			// once for the front and once for the back side of the plane
-			IncPvs(*ray->intersections[0].mObject, fpvs, bpvs, 
-				   cf, frontId, backId, frontAndBackId);
-		}
-		// the source object in the origin of the ray
-		if (ray->sourceObject.mObject)
-		{
-			IncPvs(*ray->sourceObject.mObject, fpvs, bpvs, 
-				   cf, frontId, backId, frontAndBackId);
-		}
-
-	}
-
+	if (0)
 	Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
 		  << " frontpvs: " << frontPvs << " pFront: " << pFront 
 		  << " backpvs: " << backPvs << " pBack: " << pBack 
-		  << " val " << val << " new size: " << ComputePvsSize(rays)
-		  << " fp " << fpvs << " bp " << bpvs << endl << endl;
-
-}
+		  << " val " << val << " new size: " << ComputePvsSize(rays)<< endl << endl;
+
 	return val;
 }
@@ -2016,14 +2039,20 @@
 		BoundedRay *bRay = rays.back();
 		Ray *ray = bRay->mRay;
+		float minT = bRay->mMinT;
+		float maxT = bRay->mMaxT;
+
 		rays.pop_back();
+	
+		Vector3 entP, extP;
+
+		const int cf = 
+			ray->ClassifyPlane(plane, minT, maxT, entP, extP);
 		
-		const int cf = 
-			ray->ClassifyPlane(plane, bRay->mMinT, bRay->mMaxT);
-		
+		// set id to ray classification
 		ray->SetId(cf);
 
 		switch (cf)
 		{
-		case Ray::COINCIDENT:
+		case Ray::COINCIDENT: // TODO: should really discard ray?
 			//frontRays.push_back(bRay);
 			DEL_PTR(bRay);
@@ -2038,11 +2067,10 @@
 			{
 				// find intersection of ray segment with plane
-				const Vector3 extp = ray->Extrap(bRay->mMaxT);
-				const float t = plane.FindT(ray->GetLoc(), extp);
+				const float t = plane.FindT(ray->GetLoc(), extP);
 				
-				const float newT = t * bRay->mMaxT;
-
-				frontRays.push_back(new BoundedRay(ray, bRay->mMinT, newT));
-				backRays.push_back(new BoundedRay(ray, newT, bRay->mMaxT));
+				const float newT = t * maxT;
+
+				frontRays.push_back(new BoundedRay(ray, minT, newT));
+				backRays.push_back(new BoundedRay(ray, newT, maxT));
 
 				DEL_PTR(bRay);
@@ -2054,6 +2082,5 @@
 			{
 				// find intersection of ray segment with plane
-				const Vector3 extp = ray->Extrap(bRay->mMaxT);
-				const float t = plane.FindT(ray->GetLoc(), extp);
+				const float t = plane.FindT(ray->GetLoc(), extP);
 				const float newT = t * bRay->mMaxT;
 
@@ -2411,5 +2438,5 @@
 	for (int i = 0; i < (int)mPolys.size(); ++ i)
 	{
-		const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001);
+		const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
 			
 		// split new polygon with all previous planes
@@ -2473,5 +2500,5 @@
 
 		const int cf = 
-			planePoly->ClassifyPlane(plane, 0.00001);
+			planePoly->ClassifyPlane(plane, 0.00001f);
 			
 		// split new polygon with all previous planes
Index: trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp	(revision 406)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp	(revision 406)
@@ -0,0 +1,90 @@
+#include "RenderSimulator.h"
+#include "KdTree.h"
+#include "ViewCellBsp.h"
+#include "ViewCell.h"
+
+RenderSimulator::RenderSimulator()
+{}
+
+RenderSimulator::RenderSimulator(float objRenderCost, float vcOverhead):
+mObjRenderCost(objRenderCost), mVcOverhead(vcOverhead)
+{
+}
+
+/*****************************************************
+ *     class ViewCellRenderSimulator implementation  *
+ *****************************************************/
+
+BspTreeRenderSimulator::BspTreeRenderSimulator(float objRenderCost, 
+											   float vcOverhead, 
+											   BspTree *bspTree):
+RenderSimulator(objRenderCost, vcOverhead), mBspTree(bspTree)
+{
+}
+
+Real BspTreeRenderSimulator::SimulateRendering()
+{
+	Real renderTime = 0;
+
+	// render time for 1 object of PVS
+	const float objRt = mObjRenderCost; 
+	 // const overhead for crossing a view cell border
+	const float vcOverhead = mVcOverhead;
+	
+	// total area of view cells
+	float totalArea = 0;//= mKdTree->GetBox().SurfaceArea();
+
+	ViewCellContainer viewCells;
+
+	mBspTree->CollectViewCells(viewCells);
+
+	ViewCellContainer::const_iterator it, it_end = viewCells.end();
+	
+	for (it = viewCells.begin(); it != it_end; ++ it)
+	{
+		// surface area substitute for probability
+		PolygonContainer cell;
+		float area = 0;
+
+		mBspTree->ConstructGeometry(dynamic_cast<BspViewCell *>(*it), cell);
+
+		area = Polygon3::GetArea(cell);
+		
+		renderTime += area * RenderPvs(*(*it), objRt) + vcOverhead;
+		totalArea += area;
+		CLEAR_CONTAINER(cell);
+	}
+
+	renderTime /= totalArea;
+
+	//Debug << "render time without overhead: " << renderTime * 1e-3 << endl;
+	//renderTime += (float)viewCells.size() * vcOverhead;
+
+	return renderTime;
+}
+
+Real BspTreeRenderSimulator::RenderPvs(ViewCell &viewCell, float objRenderTime) const
+{
+	return viewCell.GetPvs().GetSize() * objRenderTime;
+}
+
+/*****************************************************
+ *     class KdLeafRenderSimulator implementation    *
+ *****************************************************/
+
+KdTreeRenderSimulator::KdTreeRenderSimulator(float objRenderCost, 
+											 float vcOverhead, 
+											 KdTree *kdTree):
+RenderSimulator(objRenderCost, vcOverhead), mKdTree(kdTree)
+{
+}
+
+Real KdTreeRenderSimulator::SimulateRendering()
+{
+	Real renderTime = 0;
+
+	// total area of view cells
+	float totalArea = 0;//= mKdTree->GetBox().SurfaceArea();
+
+	return renderTime;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 369)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 406)
@@ -172,35 +172,37 @@
 int
 Ray::ClassifyPlane(const Plane3 &plane, 
-									 const float minT, 
-									 const float maxT) const
-{
-	const Vector3 entp = Extrap(minT);
-	const Vector3 extp = Extrap(maxT);
-  
-	const int entSide = plane.Side(entp);
-	const int extSide = plane.Side(extp);
+				   const float minT, 
+				   const float maxT,
+				   Vector3 &entP,
+				   Vector3 &extP) const
+{
+	entP = Extrap(minT);
+	extP = Extrap(maxT);
+  
+	const int entSide = plane.Side(entP);
+	const int extSide = plane.Side(extP);
 
 	if (entSide < 0)
+	{
+		if (extSide >= 0)
 		{
-			if (extSide >= 0)
-				{
-					return BACK_FRONT;
-				}
-			return BACK;
+			return BACK_FRONT;
 		}
+		return BACK;
+	}
 	else if (entSide > 0)
-		{
-			if (extSide <= 0)
-				return FRONT_BACK;
+	{
+		if (extSide <= 0)
+			return FRONT_BACK;
 			
-			return FRONT;
-		}
+		return FRONT;
+	}
 	else if (entSide == 0)
-		{
-			if (extSide > 0)
-				return BACK_FRONT;
-			else if (extSide < 0)
-				return FRONT_BACK;
-		}
+	{
+		if (extSide > 0)
+			return BACK_FRONT;
+		else if (extSide < 0)
+			return FRONT_BACK;
+	}
 	
 	return COINCIDENT;
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.h	(revision 406)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.h	(revision 406)
@@ -0,0 +1,513 @@
+#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/SamplingPreprocessor.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h	(revision 392)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h	(revision 406)
@@ -99,17 +99,9 @@
 	int PostprocessViewCells(const RayContainer &rays);
 
-	/** Simulated rendering using a simple render heuristics. Used to evaluate the 
-		quality of the view cell partition.
-	*/
-	Real SimulateRendering();
-	/** Simulates rendering of the pvs of one view cell, with given rendering time for an object.
-		@param viewCell the view cell holding the Pvs
-		@param objRenderTime estimated render time for one object of the Pvs
-	*/
-	Real RenderPvs(ViewCell &viewCell, const float objRenderTime) const;
-
+	// rendering costs for a single object
+	float mObjRenderCost;
+	// rendering overhead when crossing a view cell border
+	float mVcOverhead;
 };
 
-
-
 #endif
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.cpp	(revision 406)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.cpp	(revision 406)
@@ -0,0 +1,840 @@
+#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;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 404)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 406)
@@ -7,4 +7,5 @@
 #include "Polygon3.h"
 #include "ViewCell.h"
+#include "RenderSimulator.h"
 
 SamplingPreprocessor::SamplingPreprocessor(): mPass(0), mSampleRays(NULL)
@@ -16,4 +17,6 @@
   environment->GetIntValue("ViewCells.PostProcessing.samples", mPostProcessSamples);
   environment->GetIntValue("BspTree.Visualization.samples", mVisualizationSamples);
+  environment->GetFloatValue("Simulation.objRenderCost", mObjRenderCost);
+  environment->GetFloatValue("Simulation.vcOverhead", mVcOverhead);
 
   mKdPvsDepth = 100;
@@ -540,5 +543,5 @@
 		//-- render simulation
 		cout << "\nevaluating render time before merge ... ";
-		Real rt = SimulateRendering();
+		Real rt = BspTreeRenderSimulator(mObjRenderCost, mVcOverhead, mBspTree).SimulateRendering();
 		
 		cout << "avg render time: " << rt * 1e-3 << endl;
@@ -586,5 +589,5 @@
 		cout << "\nevaluating render time after merge ... ";
 			
-		rt = SimulateRendering();
+		rt = BspTreeRenderSimulator(mObjRenderCost, mVcOverhead, mBspTree).SimulateRendering();
 
 		cout << "render time: " << rt * 1e-3 << endl;
@@ -1034,50 +1037,2 @@
 	}
 }
-
-
-Real SamplingPreprocessor::RenderPvs(ViewCell &viewCell, 
-									 const float objRenderTime) const
-{
-	return viewCell.GetPvs().GetSize() * objRenderTime;
-}
-
-Real SamplingPreprocessor::SimulateRendering()
-{
-	Real renderTime = 0;
-
-	// render time for 1 object of PVS
-	const float objRt = 1.0f; 
-	 // const overhead for crossing a view cell border
-	const float vcOverhead = 0;//0.01f;
-	
-	// total area of view cells
-	float totalArea = 0;//= mKdTree->GetBox().SurfaceArea();
-
-	ViewCellContainer viewCells;
-
-	mBspTree->CollectViewCells(viewCells);
-
-	ViewCellContainer::const_iterator it, it_end = viewCells.end();
-	
-	for (it = viewCells.begin(); it != it_end; ++ it)
-	{
-		// surface area substitute for probability
-		PolygonContainer cell;
-		float area = 0;
-
-		mBspTree->ConstructGeometry(dynamic_cast<BspViewCell *>(*it), cell);
-
-		area = Polygon3::GetArea(cell);
-		
-		renderTime += area * RenderPvs(*(*it), objRt) + vcOverhead;
-		totalArea += area;
-		CLEAR_CONTAINER(cell);
-	}
-
-	renderTime /= totalArea;
-
-	//Debug << "render time without overhead: " << renderTime * 1e-3 << endl;
-	//renderTime += (float)viewCells.size() * vcOverhead;
-
-	return renderTime;
-}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsKd.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsKd.h	(revision 405)
+++ 	(revision )
@@ -1,518 +1,0 @@
-#ifndef _KdTree_H__
-#define _KdTree_H__
-
-#include <functional>
-using namespace std;
-
-#include "Containers.h"
-#include "AxisAlignedBox3.h"
-#include "Ray.h"
-#include "Pvs.h"
-
-  
-class KdNode;
-class KdLeaf;
-class KdInterior;
-class Intersectable;
-
-
-// --------------------------------------------------------------
-// Static statistics for kd-tree search
-// --------------------------------------------------------------
-class KdTreeStatistics 
-{
-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
-  KdTreeStatistics() {
-    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 KdTreeStatistics &stat) {
-    stat.Print(s);
-    return s;
-  }
-  
-};
-
-  
-class KdInterior;
-/** Abstract class for kd-tree node */
-class KdNode {
-public:
-  static int mailID;
-  int mailbox;
-  
-  void Mail() { mailbox = mailID; }
-  static void NewMail() { mailID++; }
-  bool Mailed() const { return mailbox == mailID; }
-
-
-  KdNode(KdInterior *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 */
-  KdInterior *mParent;
-  int mDepth;
-};
-
-/** Implementation of the kd-tree interior node */
-class KdInterior : public KdNode {
-
-public:
-    
-  KdInterior(KdInterior *parent):KdNode(parent), mBack(NULL), mFront(NULL) {}
-
-  /** \sa KdNode::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 */
-  KdNode *mBack;
-  /** front node */
-  KdNode *mFront;
-
-  void SetupChildLinks(KdNode *b, KdNode *f) {
-    mBack = b;
-    mFront = f;
-    b->mParent = f->mParent = this;
-  }
-  
-  void ReplaceChildLink(KdNode *oldChild, KdNode *newChild) {
-    if (mBack == oldChild)
-      mBack = newChild;
-    else
-      mFront = newChild;
-  }
-  
-  
-};
-  
-  
-/** Implementation of the kd-tree leaf node */
-class KdLeaf : public KdNode {
-public:
-  KdLeaf(KdInterior *parent, const int objects):KdNode(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 KdNode::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 KdTree nodes */
-  KdPvs mKdPvs;
-	
-	/** PVS consisting of visible objects */
-  ViewCellPvs mPvs;
-};
-
-  
-
-/** KdTree for indexing scene entities - occluders/occludees/viewcells */
-class KdTree {
-
-protected:
-  struct TraversalData
-  {  
-    KdNode *mNode;
-    AxisAlignedBox3 mBox;
-    int mDepth;
-    float mPriority;
-    
-    TraversalData() {}
-
-    TraversalData(KdNode *n, const float p):
-      mNode(n), mPriority(p)
-    {}
-    
-    TraversalData(KdNode *n,
-		   const AxisAlignedBox3 &b,
-		   const int d):
-      mNode(n), mBox(b), mDepth(d) {}
-    
-
-    bool operator<(
-		   const TraversalData &b) const {
-      KdLeaf *leafa = (KdLeaf *) mNode;
-      KdLeaf *leafb = (KdLeaf *) 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};
-  
-  KdTree();
-  
-    
-  /** 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 KdNode *Subdivide(const TraversalData &tdata);
-
-  /** Get the root of the tree */
-  KdNode *GetRoot() const {
-    return mRoot;
-  }
-
-
-  AxisAlignedBox3 GetBox() const { return mBox; }
-
-  int
-  CastRay(
-	  Ray &ray
-	  );
-
-  const KdTreeStatistics &GetStatistics() const {
-    return mStat;
-  }
-
-  void
-  CollectObjects(KdNode *n, ObjectContainer &objects);
-	
-  void
-  CollectLeaves(vector<KdLeaf *> &leaves);
-	
-  AxisAlignedBox3 GetBox(const KdNode *node) const {
-    KdInterior *parent = node->mParent;
-    if (parent == NULL)
-      return mBox;
-    
-    if (!node->IsLeaf())
-      return ((KdInterior *)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;
-  }
-	
-  KdNode *
-  FindRandomNeighbor(KdNode *n,
-										 bool onlyUnmailed
-										 );
-  
-  KdNode *
-  KdTree::GetRandomLeaf(const Plane3 &halfspace);
-
-	KdNode *
-	GetRandomLeaf(const bool onlyUnmailed = false);
-
-  int
-  FindNeighbors(KdNode *n,
-		vector<KdNode *> &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 {
-    KdNode *mNode;
-    Vector3 mExitPoint;
-    float mMaxT;
-    
-    RayTraversalData() {}
-    RayTraversalData(KdNode *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(
-		KdLeaf *node,
-		const AxisAlignedBox3 &box,
-		const int axis,
-		float &position,
-		int &objectsBack,
-		int &objectsFront
-		);
-
-  void
-  SortSplitCandidates(
-		      KdLeaf *node,
-		      const int axis
-		      );
-
-  void
-  EvaluateLeafStats(const TraversalData &data);
-
-  KdNode *
-  SubdivideNode(
-		KdLeaf *leaf,
-		const AxisAlignedBox3 &box,
-		AxisAlignedBox3 &backBBox,
-		AxisAlignedBox3 &frontBBox
-		);
-
-  bool
-  TerminationCriteriaMet(const KdLeaf *leaf);
-  
-  int
-  SelectPlane(KdLeaf *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
-  KdNode *mRoot;
-  /// bounding box of the tree root
-  AxisAlignedBox3 mBox;
-  KdTreeStatistics mStat;
-  
-};
-  
-
-
-
-
-
-#endif
