Index: GTP/trunk/Lib/Vis/Preprocessing/scripts/default.env
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/scripts/default.env	(revision 651)
+++ GTP/trunk/Lib/Vis/Preprocessing/scripts/default.env	(revision 652)
@@ -289,5 +289,5 @@
 	Construction {
 		samples 50000
-		epsilon 0.00005
+		epsilon 0.0005
 		randomize false
 		renderCostWeight 1.0
Index: GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.h	(revision 648)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.h	(revision 652)
@@ -159,5 +159,39 @@
 	
 	typedef std::priority_queue<VspBspTraversalData> VspBspTraversalQueue;
-	//typedef std::stack<VspBspTraversalData> VspBspTraversalQueue;
+	
+	
+	struct VspBspSplitCandidate
+	{  
+		/// the current node
+		Plane3 mSplitPlane;
+		
+		// parent data
+		VspBspTraversalData mParentData;
+
+		float mRenderCost;
+
+		VspBspSplitCandidate(const Plane3 &plane, const VspBspTraversalData &tData): 
+		mSplitPlane(plane), mParentData(tData)
+		{}
+
+		/** Returns cost of the traversal data.
+		*/
+		float GetCost() const
+		{
+#if 1
+			return mRenderCost;
+#endif
+#if 0
+			return (float) (-mDepth); // for kd tree
+#endif
+		}
+
+		friend bool operator<(const VspBspSplitCandidate &a, const VspBspSplitCandidate &b)
+		{
+			return a.GetCost() < b.GetCost();
+		}
+    };
+
+	typedef std::priority_queue<VspBspSplitCandidate> VspBspSplitQueue;
 
 	/** Default constructor creating an empty tree.
@@ -355,4 +389,9 @@
 								   float &pFront,
 								   float &pBack) const;
+
+	void EvalSplitCandidate(VspBspTraversalData &tData, VspBspSplitCandidate &splitData);
+
+	float EvalRenderCostDecrease(const Plane3 &candidatePlane,
+								 const VspBspTraversalData &data) const;
 
 	/** Returns view cell corresponding to 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.cpp	(revision 650)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.cpp	(revision 652)
@@ -643,4 +643,24 @@
 
 
+void VspBspTree::EvalSplitCandidate(VspBspTraversalData &tData,
+									VspBspSplitCandidate &splitData)
+{
+	VspBspTraversalData frontData;
+	VspBspTraversalData backData;
+
+	int splitAxis = 0;
+
+	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
+	SelectPlane(splitData.mSplitPlane, leaf, tData, frontData, backData, splitAxis);
+
+	// TODO: reuse
+	DEL_PTR(frontData.mGeometry);
+	DEL_PTR(backData.mGeometry);
+	
+	splitData.mRenderCost = EvalRenderCostDecrease(splitData.mSplitPlane, tData);
+	splitData.mParentData = tData;
+}
+
+
 BspNode *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
 								   VspBspTraversalData &frontData,
@@ -653,10 +673,11 @@
 	Plane3 splitPlane;
 	
-	int maxCostMisses = tData.mMaxCostMisses;
-
 	int splitAxis;
 
 	const bool success = 
 		SelectPlane(splitPlane, leaf, tData, frontData, backData, splitAxis);
+
+
+	int maxCostMisses = tData.mMaxCostMisses;
 
 	if (!success)
@@ -672,5 +693,64 @@
 	}
 
-	//! error also computed if cost ratio is missed
+	//-- the front and back traversal data is filled with the new values
+	frontData.mDepth = tData.mDepth + 1;
+	frontData.mPolygons = new PolygonContainer();
+	frontData.mRays = new RayInfoContainer();
+	
+	backData.mDepth = tData.mDepth + 1;
+	backData.mPolygons = new PolygonContainer();
+	backData.mRays = new RayInfoContainer();
+	
+	// subdivide rays
+	SplitRays(splitPlane,
+			  *tData.mRays,
+			  *frontData.mRays,
+			  *backData.mRays);
+
+
+	// how often was max cost ratio missed in this branch?
+	frontData.mMaxCostMisses = maxCostMisses;
+	backData.mMaxCostMisses = maxCostMisses;
+
+	// compute pvs
+	frontData.mPvs = ComputePvsSize(*frontData.mRays);
+	backData.mPvs = ComputePvsSize(*backData.mRays);
+
+	// split front and back node geometry and compute area
+	
+	// if geometry was not already computed
+	if (!frontData.mGeometry && !backData.mGeometry)
+	{
+		frontData.mGeometry = new BspNodeGeometry();
+		backData.mGeometry = new BspNodeGeometry();
+
+		tData.mGeometry->SplitGeometry(*frontData.mGeometry,
+									   *backData.mGeometry,
+									   splitPlane,
+									   mBox,
+									   mEpsilon);
+		
+		if (mUseAreaForPvs)
+		{
+			frontData.mProbability = frontData.mGeometry->GetArea();
+			backData.mProbability = backData.mGeometry->GetArea();
+		}
+		else
+		{
+			frontData.mProbability = frontData.mGeometry->GetVolume();
+			backData.mProbability =  tData.mProbability - frontData.mProbability;
+		}
+	}
+	
+	
+    // subdivide polygons
+	SplitPolygons(splitPlane,
+				  *tData.mPolygons,
+	              *frontData.mPolygons,
+				  *backData.mPolygons,
+				  coincident);
+
+
+
 	if (splitAxis < 3)
 		++ mBspStats.splits[splitAxis];
@@ -680,4 +760,8 @@
 	mBspStats.nodes += 2;
 
+
+
+	///////////////////////////////////////
+	
 	//-- subdivide further
 	BspInterior *interior = new BspInterior(splitPlane);
@@ -687,62 +771,4 @@
 #endif
 
-	//-- the front and back traversal data is filled with the new values
-	frontData.mDepth = tData.mDepth + 1;
-	frontData.mPolygons = new PolygonContainer();
-	frontData.mRays = new RayInfoContainer();
-	
-	backData.mDepth = tData.mDepth + 1;
-	backData.mPolygons = new PolygonContainer();
-	backData.mRays = new RayInfoContainer();
-	
-	// subdivide rays
-	SplitRays(interior->GetPlane(),
-			  *tData.mRays,
-			  *frontData.mRays,
-			  *backData.mRays);
-
-	// subdivide polygons
-	SplitPolygons(interior->GetPlane(),
-				  *tData.mPolygons,
-	              *frontData.mPolygons,
-				  *backData.mPolygons,
-				  coincident);
-
-
-	// how often was max cost ratio missed in this branch?
-	frontData.mMaxCostMisses = maxCostMisses;
-	backData.mMaxCostMisses = maxCostMisses;
-
-	// compute pvs
-	frontData.mPvs = ComputePvsSize(*frontData.mRays);
-	backData.mPvs = ComputePvsSize(*backData.mRays);
-
-	// split front and back node geometry and compute area
-	
-	// if geometry was not already computed
-	if (!frontData.mGeometry && !backData.mGeometry)
-	{
-		frontData.mGeometry = new BspNodeGeometry();
-		backData.mGeometry = new BspNodeGeometry();
-
-		tData.mGeometry->SplitGeometry(*frontData.mGeometry,
-									   *backData.mGeometry,
-									   interior->GetPlane(),
-									   mBox,
-									   mEpsilon);
-		
-		if (mUseAreaForPvs)
-		{
-			frontData.mProbability = frontData.mGeometry->GetArea();
-			backData.mProbability = backData.mGeometry->GetArea();
-		}
-		else
-		{
-			frontData.mProbability = frontData.mGeometry->GetVolume();
-			backData.mProbability =  tData.mProbability - frontData.mProbability;
-		}
-	}
-	
-
 	//-- create front and back leaf
 
@@ -767,7 +793,5 @@
 
 	interior->mTimeStamp = mTimeStamp ++;
-	//frontData.mNode->mTimeStamp = ++ mTimeStamp;
-	//backData.mNode->mTimeStamp = mTimeStamp;
-
+	
 	//DEL_PTR(leaf);
 	return interior;
@@ -1135,5 +1159,5 @@
 bool VspBspTree::SelectPlane(Plane3 &bestPlane,
 							 BspLeaf *leaf,
-							 VspBspTraversalData &data,
+							 VspBspTraversalData &data,							 
 							 VspBspTraversalData &frontData,
 							 VspBspTraversalData &backData,
@@ -1345,4 +1369,72 @@
 
 
+float VspBspTree::EvalRenderCostDecrease(const Plane3 &candidatePlane,
+										 const VspBspTraversalData &data) const
+{
+	int pvsFront = 0;
+	int pvsBack = 0;
+	int totalPvs = data.mPvs;
+
+	// probability that view point lies in back / front node
+	float pOverall = data.mProbability;
+	float pFront = 0;
+	float pBack = 0;
+
+
+	// create unique ids for pvs heuristics
+	GenerateUniqueIdsForPvs();
+	
+	for (int i = 0; i < data.mRays->size(); ++ i)
+	{
+		RayInfo rayInf = (*data.mRays)[i];
+
+		float t;
+		VssRay *ray = rayInf.mRay;
+		const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
+
+		// find front and back pvs for origing and termination object
+		AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
+		AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
+	}
+
+
+	BspNodeGeometry geomFront;
+	BspNodeGeometry geomBack;
+
+	// construct child geometry with regard to the candidate split plane
+	data.mGeometry->SplitGeometry(geomFront,
+								  geomBack,
+								  candidatePlane,
+								  mBox,
+								  mEpsilon);
+
+	if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
+	{
+		pFront = geomFront.GetVolume();
+		pBack = pOverall - pFront;
+	}
+	else
+	{
+		pFront = geomFront.GetArea();
+		pBack = geomBack.GetArea();
+	}
+	
+
+	// -- pvs rendering heuristics
+	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
+	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
+
+	// only render cost heuristics or combined with standard deviation
+	const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
+    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
+	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
+			
+	const float oldRenderCost = pOverall * penaltyOld;
+	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
+
+	return oldRenderCost - newRenderCost;
+}
+
+
 float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
 									 const VspBspTraversalData &data,
@@ -1354,10 +1446,8 @@
 	float cost = 0;
 
-	float sumBalancedRays = 0;
-	float sumRaySplits = 0;
-
+	int totalPvs = 0;
 	int pvsFront = 0;
 	int pvsBack = 0;
-
+	
 	// probability that view point lies in back / front node
 	float pOverall = 0;
@@ -1365,10 +1455,10 @@
 	pBack = 0;
 
-	int raysFront = 0;
-	int raysBack = 0;
-	int totalPvs = 0;
 
 	int limit;
 	bool useRand;
+
+	// create unique ids for pvs heuristics
+	GenerateUniqueIdsForPvs();
 
 	// choose test rays randomly if too much
@@ -1394,105 +1484,66 @@
 		const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
 
-        if (0)
-		{
-			if (cf >= 0)
-				++ raysFront;
-			if (cf <= 0)
-				++ raysBack;
-		}
-
-		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		{
-			sumBalancedRays += cf;
-		}
-
-		if (mSplitPlaneStrategy & BALANCED_RAYS)
-		{
-			if (cf == 0)
-				++ sumRaySplits;
-		}
-
-		if (mSplitPlaneStrategy & PVS)
-		{
-			// find front and back pvs for origing and termination object
-			AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
-			AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
-		}
-	}
-
-	const float raysSize = (float)data.mRays->size() + Limits::Small;
-
-	if (mSplitPlaneStrategy & PVS)
-	{
-		// create unique ids for pvs heuristics
-		GenerateUniqueIdsForPvs();
-
-		// construct child geometry with regard to the candidate split plane
-		data.mGeometry->SplitGeometry(geomFront,
-									  geomBack,
-									  candidatePlane,
-									  mBox,
-									  mEpsilon);
-
-		pOverall = data.mProbability;
-
-		if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
-		{
-			pFront = geomFront.GetVolume();
-			pBack = pOverall - pFront;
-		}
-		else
-		{
-			pFront = geomFront.GetArea();
-			pBack = geomBack.GetArea();
-		}
-	}
-
-
-	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		cost += mLeastRaySplitsFactor * sumRaySplits / raysSize;
-
-	if (mSplitPlaneStrategy & BALANCED_RAYS)
-		cost += mBalancedRaysFactor * fabs(sumBalancedRays) / raysSize;
+		// find front and back pvs for origing and termination object
+		AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
+		AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
+	}
+
+	// construct child geometry with regard to the candidate split plane
+	data.mGeometry->SplitGeometry(geomFront,
+								  geomBack,
+								  candidatePlane,
+								  mBox,
+								  mEpsilon);
+
+	pOverall = data.mProbability;
+
+	if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
+	{
+		pFront = geomFront.GetVolume();
+		pBack = pOverall - pFront;
+	}
+	else
+	{
+		pFront = geomFront.GetArea();
+		pBack = geomBack.GetArea();
+	}
+	
 
 	// -- pvs rendering heuristics
-	if (mSplitPlaneStrategy & PVS)
-	{
-		const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
-		const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
-
-		// only render cost heuristics or combined with standard deviation
-		const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
-        const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
-		const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
+	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
+	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
+
+	// only render cost heuristics or combined with standard deviation
+	const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
+    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
+	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
 			
-		const float oldRenderCost = pOverall * penaltyOld;
-		const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
-
-		float oldCost, newCost;
-
-		// only render cost
-		if (1)
-		{
-			oldCost = oldRenderCost;
-			newCost = newRenderCost;
-		}
-		else // also considering standard deviation
-		{
-			// standard deviation is difference of back and front pvs
-			const float expectedCost = 0.5f * (penaltyFront + penaltyBack);
-
-			const float newDeviation = 0.5f * 
-				fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);
-
-			const float oldDeviation = penaltyOld;
-
-			newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
-			oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
-		}
-
-		cost += mPvsFactor * newCost / (oldCost + Limits::Small);
-		
-	}
+	const float oldRenderCost = pOverall * penaltyOld;
+	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
+
+	float oldCost, newCost;
+
+	// only render cost
+	if (1)
+	{
+		oldCost = oldRenderCost;
+		newCost = newRenderCost;
+	}
+	else // also considering standard deviation
+	{
+		// standard deviation is difference of back and front pvs
+		const float expectedCost = 0.5f * (penaltyFront + penaltyBack);
+
+		const float newDeviation = 0.5f * 
+			fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);
+
+		const float oldDeviation = penaltyOld;
+
+		newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
+		oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
+	}
+
+	cost += mPvsFactor * newCost / (oldCost + Limits::Small);
+		
 
 #ifdef _DEBUG
@@ -1503,10 +1554,5 @@
 #endif
 
-	
-	// normalize cost by sum of linear factors
-	if(0)
-		return cost / (float)mCostNormalizer;
-	else
-		return cost;
+	return cost;
 }
 
