Index: trunk/VUT/work/TestCullingTerrain/TerrainFrameListener.cpp
===================================================================
--- trunk/VUT/work/TestCullingTerrain/TerrainFrameListener.cpp	(revision 418)
+++ trunk/VUT/work/TestCullingTerrain/TerrainFrameListener.cpp	(revision 421)
@@ -646,6 +646,8 @@
 {
 	// move node rather than camera so orientation is right in the visualization
-	mCamNode->yaw(mRotX);
+	mCamNode->yaw(mRotX, Ogre::Node::TS_WORLD);
+	//mCamNode->rotate(Vector3(0,1,0), mRotX, Ogre::Node::TS_WORLD);
 	mCamNode->pitch(mRotY);
+
 	mCamNode->translate(mCamNode->getLocalAxes(), mTranslateVector);
 }
@@ -1131,5 +1133,5 @@
 		toggleUseDepthPass();
 		break;
-	case KC_S:
+	case KC_H:
 		toggleShowShadows();
 		break;
@@ -1218,5 +1220,5 @@
 		mUseDemoFps = !mUseDemoFps;
 		break;
-	case KC_A:
+	case KC_Q:
 		toggleUseArbQueries();
 		break;
@@ -1312,5 +1314,5 @@
 	
 	/* Move camera forward by keypress. */
-    if (mInputDevice->isKeyDown(KC_UP))
+    if (mInputDevice->isKeyDown(KC_UP) || mInputDevice->isKeyDown(KC_W))
 	{
 		mTranslateVector.z = -mMoveScale;
@@ -1318,9 +1320,25 @@
 	}
     /* Move camera backward by keypress. */
-    if (mInputDevice->isKeyDown(KC_DOWN))
+    if (mInputDevice->isKeyDown(KC_DOWN) || mInputDevice->isKeyDown(KC_S))
     {
 		mTranslateVector.z = mMoveScale;
 		cursorPressed = true;
     }
+
+	if (mInputDevice->isKeyDown(KC_A))
+	{
+		mTranslateVector.x -= mMoveScale;
+		mTranslateVector.y += mMoveScale;
+
+		cursorPressed = true;
+	}
+
+	if (mInputDevice->isKeyDown(KC_D))
+	{
+		mTranslateVector.x += mMoveScale;
+		mTranslateVector.y -= mMoveScale;
+
+		cursorPressed = true;
+	}
 
     if (mInputDevice->isKeyDown(KC_RIGHT))
@@ -1329,5 +1347,5 @@
 		cursorPressed = true;
     }
-	
+
     if (mInputDevice->isKeyDown(KC_LEFT))
     {
Index: trunk/VUT/work/TestCullingTerrain/TerrainFrameListener.h
===================================================================
--- trunk/VUT/work/TestCullingTerrain/TerrainFrameListener.h	(revision 418)
+++ trunk/VUT/work/TestCullingTerrain/TerrainFrameListener.h	(revision 421)
@@ -405,4 +405,5 @@
 
 	bool mUseBufferedInputMouse;
+	Matrix4 rotMatrix;
 };
 
Index: trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 421)
@@ -73,5 +73,5 @@
 	
 	height 5.0
-	
+	maxViewCells 0
 	
 	PostProcessing {
Index: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp	(revision 421)
@@ -28,10 +28,4 @@
 // Static variables
 int VspKdTreeLeaf::mailID = 0;
-int VspKdTree::sTermMaxDepth = 10;
-float VspKdTree::sTermMinSize = 0.1f;
-int VspKdTree::sTermMinPvs = 10;
-int VspKdTree::sTermMinRays = 20;
-float VspKdTree::sTermMaxCostRatio = 0.1f;
-float VspKdTree::sTermMaxRayContribution = 0.1f;
 
 /// Adds object to the pvs of the front and back node
@@ -248,19 +242,16 @@
 
 
-void VspKdTree::ParseEnvironment()
-{
-	environment->GetIntValue("VspKdTree.Termination.maxDepth", sTermMaxDepth);
-	environment->GetIntValue("VspKdTree.Termination.minPvs", sTermMinPvs);
-	environment->GetIntValue("VspKdTree.Termination.minRays", sTermMinRays);
-	environment->GetFloatValue("VspKdTree.Termination.maxRayContribution", sTermMaxRayContribution);
-	environment->GetFloatValue("VspKdTree.Termination.maxCostRatio", sTermMaxCostRatio);
-
-	environment->GetFloatValue("VspKdTree.Termination.minSize", sTermMinSize);
-	sTermMinSize = sqr(sTermMinSize);
-}
-
 // Constructor
 VspKdTree::VspKdTree()
 {	  
+	environment->GetIntValue("VspKdTree.Termination.maxDepth", mTermMaxDepth);
+	environment->GetIntValue("VspKdTree.Termination.minPvs", mTermMinPvs);
+	environment->GetIntValue("VspKdTree.Termination.minRays", mTermMinRays);
+	environment->GetFloatValue("VspKdTree.Termination.maxRayContribution", mTermMaxRayContribution);
+	environment->GetFloatValue("VspKdTree.Termination.maxCostRatio", mTermMaxCostRatio);
+
+	environment->GetFloatValue("VspKdTree.Termination.minSize", mTermMinSize);
+	mTermMinSize = sqr(mTermMinSize);
+
 	environment->GetFloatValue("VspKdTree.epsilon", epsilon);
 	environment->GetFloatValue("VspKdTree.ct_div_ci", ct_div_ci);
@@ -563,5 +554,5 @@
   }
 
-	if (costRatio > sTermMaxCostRatio) 
+	if (costRatio > mTermMaxCostRatio) 
 	{
 		//		cout<<"Too big cost ratio "<<costRatio<<endl;
@@ -861,19 +852,19 @@
 	VspKdTreeLeaf *leaf = dynamic_cast<VspKdTreeLeaf *>(data.mNode);
 
-	if (data.mDepth >= sTermMaxDepth)
+	if (data.mDepth >= mTermMaxDepth)
 		++ mStat.maxDepthNodes;
   
 	//  if ( (int)(leaf->mRays.size()) < termMinCost)
 	//    stat.minCostNodes++;
-	if (leaf->GetPvsSize() < sTermMinPvs)
+	if (leaf->GetPvsSize() < mTermMinPvs)
 		++ mStat.minPvsNodes;
 
-	if (leaf->GetPvsSize() < sTermMinRays)
+	if (leaf->GetPvsSize() < mTermMinRays)
 		++ mStat.minRaysNodes;
 
-	if (0 && leaf->GetAvgRayContribution() > sTermMaxRayContribution)
+	if (0 && leaf->GetAvgRayContribution() > mTermMaxRayContribution)
 		++ mStat.maxRayContribNodes;
 	
-	if (SqrMagnitude(data.mBox.Size()) <= sTermMinSize) 
+	if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize) 
 		++ mStat.minSizeNodes;
 
@@ -882,4 +873,12 @@
 }
 
+
+inline bool VspKdTree::TerminationCriteriaMet(const VspKdTreeLeaf *leaf, 
+											  const AxisAlignedBox3 &box) const
+{
+	return ((leaf->GetPvsSize() < mTermMinPvs) || (leaf->mRays.size() < mTermMinRays) ||
+			// (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
+			 (leaf->mDepth >= mTermMaxDepth) || SqrMagnitude(box.Size()) <= mTermMinSize);
+}
 
 
@@ -889,7 +888,5 @@
 										AxisAlignedBox3 &frontBBox)
 {
-	if ( (leaf->GetPvsSize() < sTermMinPvs) || (leaf->mRays.size() < sTermMinRays) ||
-	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
-       (leaf->mDepth >= sTermMaxDepth) || SqrMagnitude(box.Size()) <= sTermMinSize ) 
+	if (TerminationCriteriaMet(leaf, box))
 	{
 
@@ -1091,5 +1088,5 @@
 	// check if we should perform a dynamic subdivision of the leaf
 	if (// leaf->mRays.size() > (unsigned)termMinCost &&
-		(leaf->GetPvsSize() >= sTermMinPvs) && 
+		(leaf->GetPvsSize() >= mTermMinPvs) && 
 		(SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
 	{
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 421)
@@ -222,10 +222,10 @@
 
 	
-	static int mailID;
-	int mailbox;
+	static int sMailId;
+	int mMailbox;
   
-	void Mail() { mailbox = mailID; }
-	static void NewMail() { mailID++; }
-	bool Mailed() const { return mailbox == mailID; }
+	void Mail() { mMailbox = sMailId; }
+	static void NewMail() { ++ sMailId; }
+	bool Mailed() const { return mMailbox == sMailId; }
 
 protected:
@@ -495,7 +495,10 @@
 	void EvaluateViewCellsStats(BspViewCellsStatistics &stat) const;
 
+	/** Parses the environment and stores the global BSP tree parameters
+	*/
+	static void ParseEnvironment();
+
 	/// BSP tree construction method
-	int mConstructionMethod;
-
+	static int sConstructionMethod;
 protected:
 
@@ -768,8 +771,4 @@
 	/// should view cells be stored or generated in the leaves?
 	bool mGenerateViewCells;
-
-		/** Parses the environment and stores the global BSP tree parameters
-	*/
-	void ParseEnvironment();
 
 	/// maximal number of polygons before subdivision termination
@@ -809,13 +808,13 @@
 
 	// factors guiding the split plane heuristics
-	float sLeastSplitsFactor;
-	float sBalancedPolysFactor;
-	float sBalancedViewCellsFactor;
-	float sVerticalSplitsFactor;
-	float sLargestPolyAreaFactor;
-	float sBlockedRaysFactor;
-	float sLeastRaySplitsFactor;
-	float sBalancedRaysFactor;
-	float sPvsFactor;
+	float mLeastSplitsFactor;
+	float mBalancedPolysFactor;
+	float mBalancedViewCellsFactor;
+	float mVerticalSplitsFactor;
+	float mLargestPolyAreaFactor;
+	float mBlockedRaysFactor;
+	float mLeastRaySplitsFactor;
+	float mBalancedRaysFactor;
+	float mPvsFactor;
 
 	//-- thresholds used for view cells merge
@@ -830,17 +829,17 @@
 		contribution for a balanced tree.
 	*/
-	static float sLeastPolySplitsTable[4];
+	static const float sLeastPolySplitsTable[4];
 	/** Evaluates split plane classification with respect to the plane's 
 		contribution for a minimum number splits in the tree.
 	*/
-	static float sBalancedPolysTable[4];
+	static const float sBalancedPolysTable[4];
 	/** Evaluates split plane classification with respect to the plane's 
 		contribution for a minimum number of ray splits.
 	*/
-	static float sLeastRaySplitsTable[5];
+	static const float sLeastRaySplitsTable[5];
 	/** Evaluates split plane classification with respect to the plane's 
 		contribution for balanced rays.
 	*/
-	static float sBalancedRaysTable[5];
+	static const float sBalancedRaysTable[5];
 
 };
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 419)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 421)
@@ -1175,4 +1175,9 @@
 		 "bspTree");
 
+  RegisterOption("ViewCells.maxViewCells",
+		 optInt,
+		 "-view_cells_max_view_cells",
+		 "0");
+
   RegisterOption("ViewCells.PostProcessing.minPvsDif",
 		 optInt,
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 421)
@@ -14,23 +14,28 @@
 #include "Plane3.h"
 
-int BspNode::mailID = 1;
+//-- static members
+
+int BspNode::sMailId = 1;
 
 /** Evaluates split plane classification with respect to the plane's 
 	contribution for a minimum number splits in the tree.
 */
-float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
+const float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
 /** Evaluates split plane classification with respect to the plane's 
 	contribution for a balanced tree.
 */
-float BspTree::sBalancedPolysTable[] = {1, -1, 0, 0};
+const float BspTree::sBalancedPolysTable[] = {1, -1, 0, 0};
 
 /** Evaluates split plane classification with respect to the plane's 
 	contribution for a minimum number of ray splits.
 */
-float BspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
+const float BspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
 /** Evaluates split plane classification with respect to the plane's 
 	contribution for balanced rays.
 */
-float BspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
+const float BspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
+
+int BspTree::sConstructionMethod = FROM_INPUT_VIEW_CELLS;
+
 
 /****************************************************************/
@@ -283,1309 +288,6 @@
 mPvsUseArea(true)
 {
-	ParseEnvironment();
 	Randomize(); // initialise random generator for heuristics
-}
- 
-const BspTreeStatistics &BspTree::GetStatistics() const 
-{
-	return mStat;
-}
-
-void BspTreeStatistics::Print(ostream &app) const
-{
-	app << "===== BspTree statistics ===============\n";
-
-	app << setprecision(4);
-
-	app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
-
-	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
-
-	app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
-
-	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
-
-	app << "#N_SPLITS ( Number of splits )\n" << splits << "\n";
-
-	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
-	maxDepthNodes * 100 / (double)Leaves() << endl;
-
-	app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
-
-	app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
-
-	app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
-
-	app << "#N_INPUT_POLYGONS (number of input polygons )\n" <<	polys << endl;
-
-	//app << "#N_PVS: " << pvs << endl;
-
-	app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
-		 (polys + splits) / (double)polys << endl;
-	
-	app << "===== END OF BspTree statistics ==========\n";
-}
-
-
-BspTree::~BspTree()
-{
-	DEL_PTR(mRoot);
-}
-
-void BspTree::InsertViewCell(ViewCell *viewCell)
-{
-	PolygonContainer *polys = new PolygonContainer();
-
-	// extract polygons that guide the split process
-	mStat.polys += AddMeshToPolygons(viewCell->GetMesh(), *polys, viewCell);
-	mBox.Include(viewCell->GetBox()); // add to BSP aabb
-
-	InsertPolygons(polys);
-}
-
-void BspTree::InsertPolygons(PolygonContainer *polys)
-{	
-	std::stack<BspTraversalData> tStack;
-	
-	// traverse existing tree or create new tree
-    if (!mRoot)
-		mRoot = new BspLeaf();
-
-	tStack.push(BspTraversalData(mRoot, polys, 0, mRootCell, new BoundedRayContainer(), 0, 
-								 mBox.SurfaceArea(), new BspNodeGeometry()));
-
-	while (!tStack.empty())
-	{
-		// filter polygons donw the tree
-		BspTraversalData tData = tStack.top();
-	    tStack.pop();
-			
-		if (!tData.mNode->IsLeaf())
-		{
-			BspInterior *interior = dynamic_cast<BspInterior *>(tData.mNode);
-
-			//-- filter view cell polygons down the tree until a leaf is reached
-			if (!tData.mPolygons->empty())
-			{
-				PolygonContainer *frontPolys = new PolygonContainer();
-				PolygonContainer *backPolys = new PolygonContainer();
-				PolygonContainer coincident;
-
-				int splits = 0;
-		
-				// split viewcell polygons with respect to split plane
-				splits += interior->SplitPolygons(*tData.mPolygons,
-												  *frontPolys, 
-												  *backPolys,
-												  coincident);
-				
-				// extract view cells associated with the split polygons
-				ViewCell *frontViewCell = mRootCell;
-				ViewCell *backViewCell = mRootCell;
-			
-				BspTraversalData frontData(interior->GetFront(), 
-										   frontPolys, 
-										   tData.mDepth + 1,
-										   mRootCell,	
-										   tData.mRays,
-										   tData.mPvs,
-										   mBox.SurfaceArea(),
-										   new BspNodeGeometry());
-
-				BspTraversalData backData(interior->GetBack(), 
-										  backPolys,
-										  tData.mDepth + 1,
-										  mRootCell,	
-										  tData.mRays,
-										  tData.mPvs,
-										  mBox.SurfaceArea(),
-										  new BspNodeGeometry());
-
-				if (!mGenerateViewCells)
-				{
-					ExtractViewCells(frontData,
-									 backData,
-									 coincident,
-									 interior->mPlane);
-				}
-
-				// don't need coincident polygons anymore
-				CLEAR_CONTAINER(coincident);
-
-				mStat.splits += splits;
-
-				// push the children on the stack
-				tStack.push(frontData);
-				tStack.push(backData);
-			}
-
-			// cleanup
-			DEL_PTR(tData.mPolygons);
-		}
-		else
-		{
-			// reached leaf => subdivide current viewcell
-			BspNode *subRoot = Subdivide(tStack, tData);
-		}
-	}
-}
-
-int BspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
-{
-	FaceContainer::const_iterator fi;
-	
-	// copy the face data to polygons
-	for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
-	{
-		Polygon3 *poly = new Polygon3((*fi), mesh);
-		
-		if (poly->Valid())
-		{
-			poly->mParent = parent; // set parent intersectable
-			polys.push_back(poly);
-		}
-		else
-			DEL_PTR(poly);
-	}
-	return (int)mesh->mFaces.size();
-}
-
-int BspTree::AddToPolygonSoup(const ViewCellContainer &viewCells, 
-							  PolygonContainer &polys, 
-							  int maxObjects)
-{
-	int limit = (maxObjects > 0) ? 
-		Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
-  
-	int polysSize = 0;
-
-	for (int i = 0; i < limit; ++ i)
-	{
-		if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
-		{
-			mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
-			polysSize += AddMeshToPolygons(viewCells[i]->GetMesh(), polys, viewCells[i]);
-		}
-	}
-
-	return polysSize;
-}
-
-int BspTree::AddToPolygonSoup(const ObjectContainer &objects, PolygonContainer &polys, int maxObjects)
-{
-	int limit = (maxObjects > 0) ? Min((int)objects.size(), maxObjects) : (int)objects.size();
-  
-	for (int i = 0; i < limit; ++i)
-	{
-		Intersectable *object = objects[i];//*it;
-		Mesh *mesh = NULL;
-
-		switch (object->Type()) // extract the meshes
-		{
-		case Intersectable::MESH_INSTANCE:
-			mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
-			break;
-		case Intersectable::VIEW_CELL:
-			mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
-			break;
-			// TODO: handle transformed mesh instances
-		default:
-			Debug << "intersectable type not supported" << endl;
-			break;
-		}
-		
-        if (mesh) // copy the mesh data to polygons
-		{
-			mBox.Include(object->GetBox()); // add to BSP tree aabb
-			AddMeshToPolygons(mesh, polys, mRootCell);
-		}
-	}
-
-	return (int)polys.size();
-}
-
-void BspTree::Construct(const ViewCellContainer &viewCells)
-{
-	mStat.nodes = 1;
-	mBox.Initialize(); 	// initialise bsp tree bounding box
-
- 	// copy view cell meshes into one big polygon soup
-	PolygonContainer *polys = new PolygonContainer();
-	mStat.polys = AddToPolygonSoup(viewCells, *polys);
-
-	// construct tree from the view cell polygons
-	Construct(polys, new BoundedRayContainer());
-}
-
-
-void BspTree::Construct(const ObjectContainer &objects)
-{
-	mStat.nodes = 1;
-	mBox.Initialize(); 	// initialise bsp tree bounding box
-	
-	PolygonContainer *polys = new PolygonContainer();
-	
-	// copy mesh instance polygons into one big polygon soup
-	mStat.polys = AddToPolygonSoup(objects, *polys);
-
-	// construct tree from polygon soup
-	Construct(polys, new BoundedRayContainer());
-}
-
-void BspTree::Construct(const RayContainer &sampleRays)
-{
-    mStat.nodes = 1;
-	mBox.Initialize(); 	// initialise BSP tree bounding box
-	
-	PolygonContainer *polys = new PolygonContainer();
-	BoundedRayContainer *rays = new BoundedRayContainer();
-
-	RayContainer::const_iterator rit, rit_end = sampleRays.end();
-
-	long startTime = GetTime();
-
-	Debug << "**** Extracting polygons from rays ****\n";
-
-	std::map<Face *, Polygon3 *> facePolyMap;
-
-	//-- extract polygons intersected by the rays
-	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
-	{
-		Ray *ray = *rit;
-	
-		// get ray-face intersection. Store polygon representing the rays together
-		// with rays intersecting the face.
-		if (!ray->intersections.empty())
-		{
-			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->intersections[0].mObject);
-			Face *face = obj->GetMesh()->mFaces[ray->intersections[0].mFace];
-
-			std::map<Face *, Polygon3 *>::iterator it = facePolyMap.find(face);
-
-			if (it != facePolyMap.end()) 
-			{
-				//store rays if needed for heuristics
-				if (mSplitPlaneStrategy & BLOCKED_RAYS)
-					(*it).second->mPiercingRays.push_back(ray);
-			} 
-			else
-			{	//store rays if needed for heuristics
-				Polygon3 *poly = new Polygon3(face, obj->GetMesh());
-				poly->mParent = obj;
-				polys->push_back(poly);
-
-				if (mSplitPlaneStrategy & BLOCKED_RAYS)
-					poly->mPiercingRays.push_back(ray);
-
-				facePolyMap[face] = poly;
-			}
-		}
-	}
-	
-	facePolyMap.clear();
-
-	// compue bounding box
-	Polygon3::IncludeInBox(*polys, mBox);
-
-	//-- store rays
-	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
-	{
-		Ray *ray = *rit;
-        ray->SetId(-1); // reset id
-
-		float minT, maxT;
-		if (BoundRay(*ray, minT, maxT))
-			rays->push_back(new BoundedRay(ray, minT, maxT));
-	}
-
-	mStat.polys = (int)polys->size();
-
-	Debug << "**** Finished polygon extraction ****" << endl;
-	Debug << (int)polys->size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
-	Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
-
-	Construct(polys, rays);
-}
-
-void BspTree::Construct(PolygonContainer *polys, BoundedRayContainer *rays)
-{
-	std::stack<BspTraversalData> tStack;
-
-	mRoot = new BspLeaf();
-
-	BspNodeGeometry *cell = new BspNodeGeometry();
-	ConstructGeometry(mRoot, *cell);
-
-	BspTraversalData tData(mRoot, polys, 0, mRootCell, rays, 
-						   ComputePvsSize(*rays), cell->GetArea(), cell);
-
-	tStack.push(tData);
-
-	mStat.Start();
-	cout << "**** Contructing bsp tree ****\n";
-
-	while (!tStack.empty()) 
-	{
-		tData = tStack.top();
-	    tStack.pop();
-
-		// subdivide leaf node
-		BspNode *subRoot = Subdivide(tStack, tData);
-	}
-
-	mStat.Stop();
-}
-
-bool BspTree::TerminationCriteriaMet(const BspTraversalData &data) const
-{
-	return 
-		(((int)data.mPolygons->size() <= mTermMaxPolygons) || 
-		 ((int)data.mRays->size() <= mTermMaxRays) ||
-		 (data.mPvs <= mTermMinPvs) ||
-		 (data.mArea <= mTermMinArea) ||
-		 (data.mDepth >= mTermMaxDepth) ||
-		 (((float)data.mPvs / ((float)data.mRays->size() + Limits::Small)) < 
-			mTermMaxRayContribution));
-}
-
-BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
-{
-	//-- terminate traversal  
-	if (TerminationCriteriaMet(tData))		
-	{
-		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
-	
-		BspViewCell *viewCell;
-
-		// generate new view cell for each leaf
-		if (mGenerateViewCells)
-		{
-			viewCell = dynamic_cast<BspViewCell *>(ViewCell::Generate());
-		}
-		else
-		{
-			// add view cell to leaf
-			viewCell = dynamic_cast<BspViewCell *>(tData.mViewCell);
-		}
-
-		leaf->SetViewCell(viewCell);
-		viewCell->mBspLeaves.push_back(leaf);
-
-		//-- add pvs
-		if (viewCell != mRootCell)
-		{
-			int conSamp = 0, sampCon = 0;
-			leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
-			
-			mStat.contributingSamples += conSamp;
-			mStat.sampleContributions += sampCon;
-		}
-
-		EvaluateLeafStats(tData);
-		
-		//-- clean up
-		
-		// discard polygons
-		CLEAR_CONTAINER(*tData.mPolygons);
-		// discard rays
-		CLEAR_CONTAINER(*tData.mRays);
-
-		DEL_PTR(tData.mPolygons);
-		DEL_PTR(tData.mRays);
-		DEL_PTR(tData.mGeometry);
-		
-		return leaf;
-	}
-
-	//-- continue subdivision
-	PolygonContainer coincident;
-	
-	PolygonContainer *frontPolys = new PolygonContainer();
-	PolygonContainer *backPolys = new PolygonContainer();
-
-	BoundedRayContainer *frontRays = new BoundedRayContainer();
-	BoundedRayContainer *backRays = new BoundedRayContainer();
-	
-	BspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
-								new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
-	BspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
-							   new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
-
-	// create new interior node and two leaf nodes
-	BspInterior *interior = SubdivideNode(tData,
-										  tFrontData,
-										  tBackData,
-										  coincident);
-
-#ifdef _DEBUG	
-	if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
-	{	for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
-			Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
-		Debug << endl;}
-#endif
-
-	// extract view cells from coincident polygons according to plane normal
-    // only if front or back polygons are empty
-	if (!mGenerateViewCells)
-	{
-   		ExtractViewCells(tFrontData,
-						 tBackData,
-						 coincident,
-						 interior->mPlane);
-						 
-	}
-
-	// don't need coincident polygons anymory
-	CLEAR_CONTAINER(coincident);
-
-	// push the children on the stack
-	tStack.push(tFrontData);
-	tStack.push(tBackData);
-
-	// cleanup
-	DEL_PTR(tData.mNode);
-	DEL_PTR(tData.mPolygons);
-	DEL_PTR(tData.mRays);
-	DEL_PTR(tData.mGeometry);		
-
-	return interior;
-}
-
-void BspTree::ExtractViewCells(BspTraversalData &frontData,
-							   BspTraversalData &backData, 
-							   const PolygonContainer &coincident,
-							   const Plane3 splitPlane) const
-{
-	// if not empty, tree is further subdivided => don't have to find view cell
-	bool foundFront = !frontData.mPolygons->empty();
-	bool foundBack = !frontData.mPolygons->empty();
-
-	PolygonContainer::const_iterator it = 
-		coincident.begin(), it_end = coincident.end();
-
-	//-- find first view cells in front and back leafs
-	for (; !(foundFront && foundBack) && (it != it_end); ++ it)
-	{
-		if (DotProd((*it)->GetNormal(), splitPlane.mNormal) > 0)
-		{
-			backData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
-			foundBack = true;
-		}
-		else
-		{
-			frontData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
-			foundFront = true;
-		}
-	}
-}
-
-BspInterior *BspTree::SubdivideNode(BspTraversalData &tData,
-									BspTraversalData &frontData,
-									BspTraversalData &backData,
-									PolygonContainer &coincident)
-{
-	mStat.nodes += 2;
-	
-	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
-	// select subdivision plane
-	BspInterior *interior = 
-		new BspInterior(SelectPlane(leaf, tData)); 
-
-#ifdef _DEBUG
-	Debug << interior << endl;
-#endif
-	
-	// subdivide rays into front and back rays
-	SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);
-	
-	// subdivide polygons with plane
-	mStat.splits += interior->SplitPolygons(*tData.mPolygons, 
-		                                    *frontData.mPolygons, 
-											*backData.mPolygons, 
-											coincident);
-
-    // compute pvs
-	frontData.mPvs = ComputePvsSize(*frontData.mRays);
-	backData.mPvs = ComputePvsSize(*backData.mRays);
-
-	// split geometry and compute area
-	if (1)
-	{
-		tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
-									   *backData.mGeometry, 
-									   *this, 
-									   interior->mPlane);
-	
-		
-		frontData.mArea = frontData.mGeometry->GetArea();
-		backData.mArea = backData.mGeometry->GetArea();
-	}
-
-	// compute accumulated ray length
-	//frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
-	//backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
-
-	//-- create front and back leaf
-
-	BspInterior *parent = leaf->GetParent();
-
-	// replace a link from node's parent
-	if (!leaf->IsRoot())
-	{
-		parent->ReplaceChildLink(leaf, interior);
-		interior->SetParent(parent);
-	}
-	else // new root
-	{
-		mRoot = interior;
-	}
-
-	// and setup child links
-	interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
-	
-	frontData.mNode = interior->mFront;
-	backData.mNode = interior->mBack;
-	
-	//DEL_PTR(leaf);
-	return interior;
-}
-
-void BspTree::SortSplitCandidates(const PolygonContainer &polys, 
-								  const int axis, 
-								  vector<SortableEntry> &splitCandidates) const
-{
-  splitCandidates.clear();
-  
-  int requestedSize = 2 * (int)polys.size();
-  // creates a sorted split candidates array  
-  splitCandidates.reserve(requestedSize);
-  
-  PolygonContainer::const_iterator it, it_end = polys.end();
-
-  AxisAlignedBox3 box;
-
-  // insert all queries 
-  for(it = polys.begin(); it != it_end; ++ it)
-  {
-	  box.Initialize();
-	  box.Include(*(*it));
-	  
-	  splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
-      splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
-  }
-  
-  stable_sort(splitCandidates.begin(), splitCandidates.end());
-}
-
-
-float BspTree::BestCostRatio(const PolygonContainer &polys,
-							 const AxisAlignedBox3 &box,
-							 const int axis,
-							 float &position,
-							 int &objectsBack,
-							 int &objectsFront) const
-{
-	vector<SortableEntry> splitCandidates;
-
-	SortSplitCandidates(polys, axis, splitCandidates);
-	
-	// 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
-	
-	int objectsLeft = 0, objectsRight = (int)polys.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 = 1e20f;
-	vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
-
-	for(ci = splitCandidates.begin(); ci != ci_end; ++ ci) 
-	{
-		switch ((*ci).type) 
-		{
-			case SortableEntry::POLY_MIN:
-				++ objectsLeft;
-				break;
-			case SortableEntry::POLY_MAX:
-			    -- objectsRight;
-				break;
-			default:
-				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 = objectsLeft * lbox.SurfaceArea() + 
-						objectsRight * rbox.SurfaceArea();
-      
-			if (sum < minSum) 
-			{
-				minSum = sum;
-				position = (*ci).value;
-
-				objectsBack = objectsLeft;
-				objectsFront = objectsRight;
-			}
-		}
-	}
-  
-	float oldCost = (float)polys.size();
-	float newCost = mCt_div_ci + minSum / boxArea;
-	float ratio = newCost / oldCost;
-
-
-#if 0
-  Debug << "====================" << endl;
-  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
-      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
-#endif
-  return ratio;
-}
-
-bool BspTree::SelectAxisAlignedPlane(Plane3 &plane, 
-									 const PolygonContainer &polys) const
-{
-	AxisAlignedBox3 box;
-	box.Initialize();
-	
-	// create bounding box of region
-	Polygon3::IncludeInBox(polys, box);
-        
-	int objectsBack = 0, objectsFront = 0;
-	int axis = 0;
-	float costRatio = MAX_FLOAT;
-	Vector3 position;
-
-	//-- area subdivision
-	for (int i = 0; i < 3; ++ i) 
-	{
-		float p = 0;
-		float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
-		
-		if (r < costRatio)
-		{
-			costRatio = r;
-			axis = i;
-			position = p;
-		}
-	}
-	
-	if (costRatio >= mMaxCostRatio)
-		return false;
-
-	Vector3 norm(0,0,0); norm[axis] = 1.0f;
-	plane = Plane3(norm, position);
-
-	return true;
-}
-
-Plane3 BspTree::SelectPlane(BspLeaf *leaf, BspTraversalData &data)
-{
-	if (data.mPolygons->empty() && data.mRays->empty())
-	{
-		Debug << "Warning: No autopartition polygon candidate available\n";
-	
-		// return axis aligned split
-		AxisAlignedBox3 box;
-		box.Initialize();
-	
-		// create bounding box of region
-		Polygon3::IncludeInBox(*data.mPolygons, box);
-
-		const int axis = box.Size().DrivingAxis();
-		const Vector3 position = (box.Min()[axis] + box.Max()[axis])*0.5f;
-
-		Vector3 norm(0,0,0); norm[axis] = 1.0f;
-		return Plane3(norm, position);
-	}
-	
-	if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
-		((int)data.mPolygons->size() > mTermMaxPolysForAxisAligned) &&
-		((int)data.mRays->size() > mTermMaxRaysForAxisAligned) &&
-		((mTermMaxObjectsForAxisAligned < 0) || 
-		  (Polygon3::ParentObjectsSize(*data.mPolygons) > mTermMaxObjectsForAxisAligned)))
-	{
-		Plane3 plane;
-		if (SelectAxisAlignedPlane(plane, *data.mPolygons))
-			return plane;
-	}
-
-	// simplest strategy: just take next polygon
-	if (mSplitPlaneStrategy & RANDOM_POLYGON)
-	{
-        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();
-	}
-
-	// use heuristics to find appropriate plane
-	return SelectPlaneHeuristics(leaf, data); 
-}
-
-Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf, BspTraversalData &data)
-{
-	float lowestCost = MAX_FLOAT;
-	Plane3 bestPlane;
-	Plane3 plane;
-
-	int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
-	
-	int candidateIdx = limit;
-
-	for (int i = 0; i < limit; ++ i)
-	{
-		candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
-		
-		Polygon3 *poly = (*data.mPolygons)[candidateIdx];
-		// evaluate current candidate
-		const float candidateCost = 
-			SplitPlaneCost(poly->GetSupportingPlane(), data);
-
-		if (candidateCost < lowestCost)
-		{
-			bestPlane = poly->GetSupportingPlane();
-			lowestCost = candidateCost;
-		}
-	}
-	
-	//Debug << "lowest: " << lowestCost << endl;
-
-	//-- choose candidate planes extracted from rays
-	// we currently use two methods
-	// 1) take 3 ray endpoints, where two are minimum and one a maximum
-	//    point or the other way round
-	// 2) take plane normal as plane normal and the midpoint of the ray.
-	//    PROBLEM: does not resemble any point where visibility is likely to change
-	const BoundedRayContainer *rays = data.mRays;
-
-	for (int i = 0; i < mMaxRayCandidates / 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(normal, pt);
-	
-		const float candidateCost = SplitPlaneCost(plane, data);
-
-		if (candidateCost < lowestCost)
-		{
-			bestPlane = plane;
-			
-			lowestCost = candidateCost;
-		}
-	}
-
-	//Debug << "lowest: " << lowestCost << endl;
-	for (int i = 0; i < mMaxRayCandidates / 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;
-		}
-	}	
-
-#ifdef _DEBUG
-	Debug << "plane lowest cost: " << lowestCost << endl;
-#endif
-	return bestPlane;
-}
-
-int BspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
-{
-	const int candidateIdx = Random(currentIdx --);
-
-	// swap candidates to avoid testing same plane 2 times
-	std::swap(polys[currentIdx], polys[candidateIdx]);
-	
-	return currentIdx;
-	//return Random((int)polys.size());
-}
-
-float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
-							  const PolygonContainer &polys) const
-{
-	float val = 0;
-
-	float sumBalancedPolys = 0;
-	float sumSplits = 0;
-	float sumPolyArea = 0;
-	float sumBalancedViewCells = 0;
-	float sumBlockedRays = 0;
-	float totalBlockedRays = 0;
-	//float totalArea = 0;
-	int totalViewCells = 0;
-
-	// need three unique ids for each type of view cell
-	// for balanced view cells criterium
-	ViewCell::NewMail();
-	const int backId = ViewCell::sMailId;
-	ViewCell::NewMail();
-	const int frontId = ViewCell::sMailId;
-	ViewCell::NewMail();
-	const int frontAndBackId = ViewCell::sMailId;
-
-	PolygonContainer::const_iterator it, it_end = polys.end();
-
-	for (it = polys.begin(); it != it_end; ++ it)
-	{
-		const int classification = (*it)->ClassifyPlane(candidatePlane);
-
-		if (mSplitPlaneStrategy & BALANCED_POLYS)
-			sumBalancedPolys += sBalancedPolysTable[classification];
-		
-		if (mSplitPlaneStrategy & LEAST_SPLITS)
-			sumSplits += sLeastPolySplitsTable[classification];
-
-		if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
-		{
-			if (classification == Polygon3::COINCIDENT)
-				sumPolyArea += (*it)->GetArea();
-			//totalArea += area;
-		}
-		
-		if (mSplitPlaneStrategy & BLOCKED_RAYS)
-		{
-			const float blockedRays = (float)(*it)->mPiercingRays.size();
-		
-			if (classification == Polygon3::COINCIDENT)
-				sumBlockedRays += blockedRays;
-			
-			totalBlockedRays += blockedRays;
-		}
-
-		// assign view cells to back or front according to classificaion
-		if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
-		{
-			MeshInstance *viewCell = (*it)->mParent;
-		
-			// assure that we only count a view cell 
-			// once for the front and once for the back side of the plane
-			if (classification == Polygon3::FRONT_SIDE)
-			{
-				if ((viewCell->mMailbox != frontId) && 
-					(viewCell->mMailbox != frontAndBackId))
-				{
-					sumBalancedViewCells += 1.0;
-
-					if (viewCell->mMailbox != backId)
-						viewCell->mMailbox = frontId;
-					else
-						viewCell->mMailbox = frontAndBackId;
-					
-					++ totalViewCells;
-				}
-			}
-			else if (classification == Polygon3::BACK_SIDE)
-			{
-				if ((viewCell->mMailbox != backId) &&
-				    (viewCell->mMailbox != frontAndBackId))
-				{
-					sumBalancedViewCells -= 1.0;
-
-					if (viewCell->mMailbox != frontId)
-						viewCell->mMailbox = backId;
-					else
-						viewCell->mMailbox = frontAndBackId; 
-
-					++ totalViewCells;
-				}
-			}
-		}
-	}
-
-	const float polysSize = (float)polys.size() + Limits::Small;
-
-	// all values should be approx. between 0 and 1 so they can be combined
-	// and scaled with the factors according to their importance
-	if (mSplitPlaneStrategy & BALANCED_POLYS)
-		val += sBalancedPolysFactor * fabs(sumBalancedPolys) / polysSize;
-	
-	if (mSplitPlaneStrategy & LEAST_SPLITS)  
-		val += sLeastSplitsFactor * sumSplits / polysSize;
-
-	if (mSplitPlaneStrategy & LARGEST_POLY_AREA) 
-		// HACK: polys.size should be total area so scaling is between 0 and 1
-		val += sLargestPolyAreaFactor * (float)polys.size() / sumPolyArea;
-
-	if (mSplitPlaneStrategy & BLOCKED_RAYS)
-		if (totalBlockedRays != 0)
-			val += sBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;
-
-	if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
-		val += sBalancedViewCellsFactor * fabs(sumBalancedViewCells) / 
-			((float)totalViewCells + Limits::Small);
-	
-	return val;
-}
-
-bool BspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
-{
-	maxT = 1e6;
-	minT = 0;
-	
-	// test with tree bounding box
-	if (!mBox.GetMinMaxT(ray, &minT, &maxT))
-		return false;
-
-	if (minT < 0) // start ray from origin
-		minT = 0;
-
-	// bound ray or line segment
-	if ((ray.GetType() == Ray::LOCAL_RAY) && 
-	    !ray.intersections.empty() && 
-		(ray.intersections[0].mT <= maxT))
-	{
-		maxT = ray.intersections[0].mT;
-	}
-
-	return true;
-}
-
-float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
-							  const BoundedRayContainer &rays,
-							  const int pvs,
-							  const float area,
-							  const BspNodeGeometry &cell) const
-{
-	float val = 0;
-
-	float sumBalancedRays = 0;
-	float sumRaySplits = 0;
-	
-	int backId = 0;
-	int frontId = 0;
-	int frontAndBackId = 0;
-
-	int frontPvs = 0;
-	int backPvs = 0;
-
-	// probability that view point lies in child
-	float pOverall = 0;
-	float pFront = 0;
-	float pBack = 0;
-
-	if (mSplitPlaneStrategy & PVS)
-	{
-		// create three unique ids for pvs heuristics
-		Intersectable::NewMail(); backId = ViewCell::sMailId;
-		Intersectable::NewMail(); frontId = ViewCell::sMailId;
-		Intersectable::NewMail(); frontAndBackId = ViewCell::sMailId;
-
-		if (mPvsUseArea) // use front and back cell areas to approximate volume
-		{	
-			// construct child geometry with regard to the candidate split plane
-			BspNodeGeometry frontCell;
-			BspNodeGeometry backCell;
-
-			cell.SplitGeometry(frontCell, backCell, *this, candidatePlane);
-		
-			pFront = frontCell.GetArea();
-			pBack = backCell.GetArea();
-
-			pOverall = area;
-		}
-	}
-			
-	BoundedRayContainer::const_iterator rit, rit_end = rays.end();
-
-	for (rit = rays.begin(); rit != rays.end(); ++ rit)
-	{
-		Ray *ray = (*rit)->mRay;
-		const float minT = (*rit)->mMinT;
-		const float maxT = (*rit)->mMaxT;
-
-		Vector3 entP, extP;
-
-		const int cf = 
-			ray->ClassifyPlane(candidatePlane, minT, maxT, entP, extP);
-
-		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		{
-			sumBalancedRays += sBalancedRaysTable[cf];
-		}
-		
-		if (mSplitPlaneStrategy & BALANCED_RAYS)
-		{
-			sumRaySplits += sLeastRaySplitsTable[cf];
-		}
-
-		if (mSplitPlaneStrategy & PVS)
-		{
-			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, frontPvs, backPvs, 
-					   cf, frontId, backId, frontAndBackId);
-			}
-
-			// the source object in the origin of the ray
-			if (ray->sourceObject.mObject)
-			{
-				IncPvs(*ray->sourceObject.mObject, frontPvs, backPvs, 
-					   cf, frontId, backId, frontAndBackId);
-			}
-
-			if (!mPvsUseArea) // use front and back cell areas to approximate volume
-			{
-				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;
-				}
-			}
-		}
-	}
-
-	const float raysSize = (float)rays.size() + Limits::Small;
-	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		val += sLeastRaySplitsFactor * sumRaySplits / raysSize;
-
-	if (mSplitPlaneStrategy & BALANCED_RAYS)
-		val += sBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
-
-	float denom = pOverall * (float)pvs * 2.0f + Limits::Small;
-	if ((mSplitPlaneStrategy & PVS) && area && pvs)
-	{
-		val += sPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
-
-		// give penalty to unbalanced split
-		if (0)
-		if (((pFront * 0.2 + Limits::Small) > pBack) || (pFront < (pBack * 0.2 + Limits::Small)))
-			val += 0.5;
-	}
-
-#ifdef _DEBUG
-	Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
-		  << " frontpvs: " << frontPvs << " pFront: " << pFront 
-		  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
-#endif
-	return val;
-}
-
-void BspTree::IncPvs(Intersectable &obj,
-					 int &frontPvs,
-					 int &backPvs,
-					 const int cf, 
-					 const int frontId, 
-					 const int backId, 
-					 const int frontAndBackId) const
-{
-	// TODO: does this really belong to no pvs?
-	//if (cf == Ray::COINCIDENT) return;
-
-	if (cf == Ray::FRONT)
-	{
-		if ((obj.mMailbox != frontId) && 
-			(obj.mMailbox != frontAndBackId))
-		{
-			++ frontPvs;
-
-			if (obj.mMailbox != backId)
-				obj.mMailbox = frontId;
-			else
-				obj.mMailbox = frontAndBackId;					
-		}
-	}
-	else if (cf == Ray::BACK)
-	{
-		if ((obj.mMailbox != backId) &&
-			(obj.mMailbox != frontAndBackId))
-		{
-			++ backPvs;
-
-			if (obj.mMailbox != frontId)
-				obj.mMailbox = backId;
-			else
-				obj.mMailbox = frontAndBackId; 
-		}
-	}
-	// object belongs to both PVS
-	else if ((cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT) ||(cf == Ray::COINCIDENT))
-	{
-		if (obj.mMailbox !=  frontAndBackId)
-		{
-			if (obj.mMailbox != frontId)
-				++ frontPvs;
-			if (obj.mMailbox != backId)
-				++ backPvs;
-		
-			obj.mMailbox = frontAndBackId;
-		}
-	}
-}
-
-float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
-							  BspTraversalData &data) const
-{
-	float val = 0;
-
-	if (mSplitPlaneStrategy & VERTICAL_AXIS)
-	{
-		Vector3 tinyAxis(0,0,0); tinyAxis[mBox.Size().TinyAxis()] = 1.0f;
-		// we put a penalty on the dot product between the "tiny" vertical axis
-		// and the split plane axis
-		val += sVerticalSplitsFactor * 
-			   fabs(DotProd(candidatePlane.mNormal, tinyAxis));
-	}
-
-	// the following criteria loop over all polygons to find the cost value
-	if ((mSplitPlaneStrategy & BALANCED_POLYS)      ||
-		(mSplitPlaneStrategy & LEAST_SPLITS)        ||
-		(mSplitPlaneStrategy & LARGEST_POLY_AREA)   ||
-		(mSplitPlaneStrategy & BALANCED_VIEW_CELLS) ||
-		(mSplitPlaneStrategy & BLOCKED_RAYS))
-	{
-		val += SplitPlaneCost(candidatePlane, *data.mPolygons);
-	}
-
-	// the following criteria loop over all rays to find the cost value
-	if ((mSplitPlaneStrategy & BALANCED_RAYS)      ||
-		(mSplitPlaneStrategy & LEAST_RAY_SPLITS)   ||
-		(mSplitPlaneStrategy & PVS))
-	{
-		val += SplitPlaneCost(candidatePlane, *data.mRays, data.mPvs, 
-							  data.mArea, *data.mGeometry);
-	}
-
-	// return linear combination of the sums
-	return val;
-}
-
-void BspTree::ParseEnvironment()
-{
-	//-- parse bsp cell tree construction method
-	/*char constructionMethodStr[60];
-	
-	environment->GetStringValue("BspTree.Construction.input", constructionMethodStr);
-
-	mConstructionMethod = FROM_INPUT_VIEW_CELLS;
-	
-	if (strcmp(constructionMethodStr, "fromViewCells") == 0)
-		mConstructionMethod = FROM_INPUT_VIEW_CELLS;
-	else if (strcmp(constructionMethodStr, "fromSceneGeometry") == 0)
-		mConstructionMethod = FROM_SCENE_GEOMETRY;
-	else if (strcmp(constructionMethodStr, "fromRays") == 0)
-		mConstructionMethod = FROM_RAYS;
-	else 
-	{
-		cerr << "Wrong construction method " << constructionMethodStr << endl;
-		exit(1);
-    }
-
-	Debug << "Construction method: " << constructionMethodStr << endl;
-*/
+
 	//-- termination criteria for autopartition
 	environment->GetIntValue("BspTree.Termination.maxDepth", mTermMaxDepth);
@@ -1651,5 +353,1308 @@
 	if (mSplitPlaneStrategy & PVS)
 		Debug << "pvs";
+
 	Debug << endl;
+}
+
+void BspTree::ParseEnvironment()
+{
+	//-- parse bsp cell tree construction method
+	char constructionMethodStr[60];
+	
+	environment->GetStringValue("BspTree.Construction.input", constructionMethodStr);
+
+	sConstructionMethod = FROM_INPUT_VIEW_CELLS;
+	
+	if (strcmp(constructionMethodStr, "fromViewCells") == 0)
+		sConstructionMethod = FROM_INPUT_VIEW_CELLS;
+	else if (strcmp(constructionMethodStr, "fromSceneGeometry") == 0)
+		sConstructionMethod = FROM_SCENE_GEOMETRY;
+	else if (strcmp(constructionMethodStr, "fromRays") == 0)
+		sConstructionMethod = FROM_RAYS;
+	else 
+	{
+		cerr << "Wrong construction method " << constructionMethodStr << endl;
+		exit(1);
+    }
+
+	Debug << "Construction method: " << constructionMethodStr << endl;
+}
+
+const BspTreeStatistics &BspTree::GetStatistics() const 
+{
+	return mStat;
+}
+
+void BspTreeStatistics::Print(ostream &app) const
+{
+	app << "===== BspTree statistics ===============\n";
+
+	app << setprecision(4);
+
+	app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
+
+	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
+
+	app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
+
+	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
+
+	app << "#N_SPLITS ( Number of splits )\n" << splits << "\n";
+
+	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
+	maxDepthNodes * 100 / (double)Leaves() << endl;
+
+	app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
+
+	app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
+
+	app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
+
+	app << "#N_INPUT_POLYGONS (number of input polygons )\n" <<	polys << endl;
+
+	//app << "#N_PVS: " << pvs << endl;
+
+	app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
+		 (polys + splits) / (double)polys << endl;
+	
+	app << "===== END OF BspTree statistics ==========\n";
+}
+
+
+BspTree::~BspTree()
+{
+	DEL_PTR(mRoot);
+}
+
+void BspTree::InsertViewCell(ViewCell *viewCell)
+{
+	PolygonContainer *polys = new PolygonContainer();
+
+	// extract polygons that guide the split process
+	mStat.polys += AddMeshToPolygons(viewCell->GetMesh(), *polys, viewCell);
+	mBox.Include(viewCell->GetBox()); // add to BSP aabb
+
+	InsertPolygons(polys);
+}
+
+void BspTree::InsertPolygons(PolygonContainer *polys)
+{	
+	std::stack<BspTraversalData> tStack;
+	
+	// traverse existing tree or create new tree
+    if (!mRoot)
+		mRoot = new BspLeaf();
+
+	tStack.push(BspTraversalData(mRoot, polys, 0, mRootCell, new BoundedRayContainer(), 0, 
+								 mBox.SurfaceArea(), new BspNodeGeometry()));
+
+	while (!tStack.empty())
+	{
+		// filter polygons donw the tree
+		BspTraversalData tData = tStack.top();
+	    tStack.pop();
+			
+		if (!tData.mNode->IsLeaf())
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(tData.mNode);
+
+			//-- filter view cell polygons down the tree until a leaf is reached
+			if (!tData.mPolygons->empty())
+			{
+				PolygonContainer *frontPolys = new PolygonContainer();
+				PolygonContainer *backPolys = new PolygonContainer();
+				PolygonContainer coincident;
+
+				int splits = 0;
+		
+				// split viewcell polygons with respect to split plane
+				splits += interior->SplitPolygons(*tData.mPolygons,
+												  *frontPolys, 
+												  *backPolys,
+												  coincident);
+				
+				// extract view cells associated with the split polygons
+				ViewCell *frontViewCell = mRootCell;
+				ViewCell *backViewCell = mRootCell;
+			
+				BspTraversalData frontData(interior->GetFront(), 
+										   frontPolys, 
+										   tData.mDepth + 1,
+										   mRootCell,	
+										   tData.mRays,
+										   tData.mPvs,
+										   mBox.SurfaceArea(),
+										   new BspNodeGeometry());
+
+				BspTraversalData backData(interior->GetBack(), 
+										  backPolys,
+										  tData.mDepth + 1,
+										  mRootCell,	
+										  tData.mRays,
+										  tData.mPvs,
+										  mBox.SurfaceArea(),
+										  new BspNodeGeometry());
+
+				if (!mGenerateViewCells)
+				{
+					ExtractViewCells(frontData,
+									 backData,
+									 coincident,
+									 interior->mPlane);
+				}
+
+				// don't need coincident polygons anymore
+				CLEAR_CONTAINER(coincident);
+
+				mStat.splits += splits;
+
+				// push the children on the stack
+				tStack.push(frontData);
+				tStack.push(backData);
+			}
+
+			// cleanup
+			DEL_PTR(tData.mPolygons);
+		}
+		else
+		{
+			// reached leaf => subdivide current viewcell
+			BspNode *subRoot = Subdivide(tStack, tData);
+		}
+	}
+}
+
+int BspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
+{
+	FaceContainer::const_iterator fi;
+	
+	// copy the face data to polygons
+	for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
+	{
+		Polygon3 *poly = new Polygon3((*fi), mesh);
+		
+		if (poly->Valid())
+		{
+			poly->mParent = parent; // set parent intersectable
+			polys.push_back(poly);
+		}
+		else
+			DEL_PTR(poly);
+	}
+	return (int)mesh->mFaces.size();
+}
+
+int BspTree::AddToPolygonSoup(const ViewCellContainer &viewCells, 
+							  PolygonContainer &polys, 
+							  int maxObjects)
+{
+	int limit = (maxObjects > 0) ? 
+		Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
+  
+	int polysSize = 0;
+
+	for (int i = 0; i < limit; ++ i)
+	{
+		if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
+		{
+			mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
+			polysSize += AddMeshToPolygons(viewCells[i]->GetMesh(), polys, viewCells[i]);
+		}
+	}
+
+	return polysSize;
+}
+
+int BspTree::AddToPolygonSoup(const ObjectContainer &objects, PolygonContainer &polys, int maxObjects)
+{
+	int limit = (maxObjects > 0) ? Min((int)objects.size(), maxObjects) : (int)objects.size();
+  
+	for (int i = 0; i < limit; ++i)
+	{
+		Intersectable *object = objects[i];//*it;
+		Mesh *mesh = NULL;
+
+		switch (object->Type()) // extract the meshes
+		{
+		case Intersectable::MESH_INSTANCE:
+			mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
+			break;
+		case Intersectable::VIEW_CELL:
+			mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
+			break;
+			// TODO: handle transformed mesh instances
+		default:
+			Debug << "intersectable type not supported" << endl;
+			break;
+		}
+		
+        if (mesh) // copy the mesh data to polygons
+		{
+			mBox.Include(object->GetBox()); // add to BSP tree aabb
+			AddMeshToPolygons(mesh, polys, mRootCell);
+		}
+	}
+
+	return (int)polys.size();
+}
+
+void BspTree::Construct(const ViewCellContainer &viewCells)
+{
+	mStat.nodes = 1;
+	mBox.Initialize(); 	// initialise bsp tree bounding box
+
+ 	// copy view cell meshes into one big polygon soup
+	PolygonContainer *polys = new PolygonContainer();
+	mStat.polys = AddToPolygonSoup(viewCells, *polys);
+
+	// construct tree from the view cell polygons
+	Construct(polys, new BoundedRayContainer());
+}
+
+
+void BspTree::Construct(const ObjectContainer &objects)
+{
+	mStat.nodes = 1;
+	mBox.Initialize(); 	// initialise bsp tree bounding box
+	
+	PolygonContainer *polys = new PolygonContainer();
+	
+	// copy mesh instance polygons into one big polygon soup
+	mStat.polys = AddToPolygonSoup(objects, *polys);
+
+	// construct tree from polygon soup
+	Construct(polys, new BoundedRayContainer());
+}
+
+void BspTree::Construct(const RayContainer &sampleRays)
+{
+    mStat.nodes = 1;
+	mBox.Initialize(); 	// initialise BSP tree bounding box
+	
+	PolygonContainer *polys = new PolygonContainer();
+	BoundedRayContainer *rays = new BoundedRayContainer();
+
+	RayContainer::const_iterator rit, rit_end = sampleRays.end();
+
+	long startTime = GetTime();
+
+	Debug << "**** Extracting polygons from rays ****\n";
+
+	std::map<Face *, Polygon3 *> facePolyMap;
+
+	//-- extract polygons intersected by the rays
+	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
+	{
+		Ray *ray = *rit;
+	
+		// get ray-face intersection. Store polygon representing the rays together
+		// with rays intersecting the face.
+		if (!ray->intersections.empty())
+		{
+			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->intersections[0].mObject);
+			Face *face = obj->GetMesh()->mFaces[ray->intersections[0].mFace];
+
+			std::map<Face *, Polygon3 *>::iterator it = facePolyMap.find(face);
+
+			if (it != facePolyMap.end()) 
+			{
+				//store rays if needed for heuristics
+				if (mSplitPlaneStrategy & BLOCKED_RAYS)
+					(*it).second->mPiercingRays.push_back(ray);
+			} 
+			else
+			{	//store rays if needed for heuristics
+				Polygon3 *poly = new Polygon3(face, obj->GetMesh());
+				poly->mParent = obj;
+				polys->push_back(poly);
+
+				if (mSplitPlaneStrategy & BLOCKED_RAYS)
+					poly->mPiercingRays.push_back(ray);
+
+				facePolyMap[face] = poly;
+			}
+		}
+	}
+	
+	facePolyMap.clear();
+
+	// compue bounding box
+	Polygon3::IncludeInBox(*polys, mBox);
+
+	//-- store rays
+	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
+	{
+		Ray *ray = *rit;
+        ray->SetId(-1); // reset id
+
+		float minT, maxT;
+		if (BoundRay(*ray, minT, maxT))
+			rays->push_back(new BoundedRay(ray, minT, maxT));
+	}
+
+	mStat.polys = (int)polys->size();
+
+	Debug << "**** Finished polygon extraction ****" << endl;
+	Debug << (int)polys->size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
+	Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+
+	Construct(polys, rays);
+}
+
+void BspTree::Construct(PolygonContainer *polys, BoundedRayContainer *rays)
+{
+	std::stack<BspTraversalData> tStack;
+
+	mRoot = new BspLeaf();
+
+	BspNodeGeometry *cell = new BspNodeGeometry();
+	ConstructGeometry(mRoot, *cell);
+
+	BspTraversalData tData(mRoot, polys, 0, mRootCell, rays, 
+						   ComputePvsSize(*rays), cell->GetArea(), cell);
+
+	tStack.push(tData);
+
+	mStat.Start();
+	cout << "**** Contructing bsp tree ****\n";
+
+	while (!tStack.empty()) 
+	{
+		tData = tStack.top();
+	    tStack.pop();
+
+		// subdivide leaf node
+		BspNode *subRoot = Subdivide(tStack, tData);
+	}
+
+	mStat.Stop();
+}
+
+bool BspTree::TerminationCriteriaMet(const BspTraversalData &data) const
+{
+	return 
+		(((int)data.mPolygons->size() <= mTermMaxPolygons) || 
+		 ((int)data.mRays->size() <= mTermMaxRays) ||
+		 (data.mPvs <= mTermMinPvs) ||
+		 (data.mArea <= mTermMinArea) ||
+		 (data.mDepth >= mTermMaxDepth) ||
+		 (((float)data.mPvs / ((float)data.mRays->size() + Limits::Small)) < 
+			mTermMaxRayContribution));
+}
+
+BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
+{
+	//-- terminate traversal  
+	if (TerminationCriteriaMet(tData))		
+	{
+		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
+	
+		BspViewCell *viewCell;
+
+		// generate new view cell for each leaf
+		if (mGenerateViewCells)
+		{
+			viewCell = dynamic_cast<BspViewCell *>(ViewCell::Generate());
+		}
+		else
+		{
+			// add view cell to leaf
+			viewCell = dynamic_cast<BspViewCell *>(tData.mViewCell);
+		}
+
+		leaf->SetViewCell(viewCell);
+		viewCell->mBspLeaves.push_back(leaf);
+
+		//-- add pvs
+		if (viewCell != mRootCell)
+		{
+			int conSamp = 0, sampCon = 0;
+			leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
+			
+			mStat.contributingSamples += conSamp;
+			mStat.sampleContributions += sampCon;
+		}
+
+		EvaluateLeafStats(tData);
+		
+		//-- clean up
+		
+		// discard polygons
+		CLEAR_CONTAINER(*tData.mPolygons);
+		// discard rays
+		CLEAR_CONTAINER(*tData.mRays);
+
+		DEL_PTR(tData.mPolygons);
+		DEL_PTR(tData.mRays);
+		DEL_PTR(tData.mGeometry);
+		
+		return leaf;
+	}
+
+	//-- continue subdivision
+	PolygonContainer coincident;
+	
+	PolygonContainer *frontPolys = new PolygonContainer();
+	PolygonContainer *backPolys = new PolygonContainer();
+
+	BoundedRayContainer *frontRays = new BoundedRayContainer();
+	BoundedRayContainer *backRays = new BoundedRayContainer();
+	
+	BspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
+								new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
+	BspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
+							   new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
+
+	// create new interior node and two leaf nodes
+	BspInterior *interior = SubdivideNode(tData,
+										  tFrontData,
+										  tBackData,
+										  coincident);
+
+#ifdef _DEBUG	
+	if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
+	{	for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
+			Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
+		Debug << endl;}
+#endif
+
+	// extract view cells from coincident polygons according to plane normal
+    // only if front or back polygons are empty
+	if (!mGenerateViewCells)
+	{
+   		ExtractViewCells(tFrontData,
+						 tBackData,
+						 coincident,
+						 interior->mPlane);
+						 
+	}
+
+	// don't need coincident polygons anymory
+	CLEAR_CONTAINER(coincident);
+
+	// push the children on the stack
+	tStack.push(tFrontData);
+	tStack.push(tBackData);
+
+	// cleanup
+	DEL_PTR(tData.mNode);
+	DEL_PTR(tData.mPolygons);
+	DEL_PTR(tData.mRays);
+	DEL_PTR(tData.mGeometry);		
+
+	return interior;
+}
+
+void BspTree::ExtractViewCells(BspTraversalData &frontData,
+							   BspTraversalData &backData, 
+							   const PolygonContainer &coincident,
+							   const Plane3 splitPlane) const
+{
+	// if not empty, tree is further subdivided => don't have to find view cell
+	bool foundFront = !frontData.mPolygons->empty();
+	bool foundBack = !frontData.mPolygons->empty();
+
+	PolygonContainer::const_iterator it = 
+		coincident.begin(), it_end = coincident.end();
+
+	//-- find first view cells in front and back leafs
+	for (; !(foundFront && foundBack) && (it != it_end); ++ it)
+	{
+		if (DotProd((*it)->GetNormal(), splitPlane.mNormal) > 0)
+		{
+			backData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
+			foundBack = true;
+		}
+		else
+		{
+			frontData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
+			foundFront = true;
+		}
+	}
+}
+
+BspInterior *BspTree::SubdivideNode(BspTraversalData &tData,
+									BspTraversalData &frontData,
+									BspTraversalData &backData,
+									PolygonContainer &coincident)
+{
+	mStat.nodes += 2;
+	
+	BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
+	// select subdivision plane
+	BspInterior *interior = 
+		new BspInterior(SelectPlane(leaf, tData)); 
+
+#ifdef _DEBUG
+	Debug << interior << endl;
+#endif
+	
+	// subdivide rays into front and back rays
+	SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);
+	
+	// subdivide polygons with plane
+	mStat.splits += interior->SplitPolygons(*tData.mPolygons, 
+		                                    *frontData.mPolygons, 
+											*backData.mPolygons, 
+											coincident);
+
+    // compute pvs
+	frontData.mPvs = ComputePvsSize(*frontData.mRays);
+	backData.mPvs = ComputePvsSize(*backData.mRays);
+
+	// split geometry and compute area
+	if (1)
+	{
+		tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
+									   *backData.mGeometry, 
+									   *this, 
+									   interior->mPlane);
+	
+		
+		frontData.mArea = frontData.mGeometry->GetArea();
+		backData.mArea = backData.mGeometry->GetArea();
+	}
+
+	// compute accumulated ray length
+	//frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
+	//backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
+
+	//-- create front and back leaf
+
+	BspInterior *parent = leaf->GetParent();
+
+	// replace a link from node's parent
+	if (!leaf->IsRoot())
+	{
+		parent->ReplaceChildLink(leaf, interior);
+		interior->SetParent(parent);
+	}
+	else // new root
+	{
+		mRoot = interior;
+	}
+
+	// and setup child links
+	interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
+	
+	frontData.mNode = interior->mFront;
+	backData.mNode = interior->mBack;
+	
+	//DEL_PTR(leaf);
+	return interior;
+}
+
+void BspTree::SortSplitCandidates(const PolygonContainer &polys, 
+								  const int axis, 
+								  vector<SortableEntry> &splitCandidates) const
+{
+  splitCandidates.clear();
+  
+  int requestedSize = 2 * (int)polys.size();
+  // creates a sorted split candidates array  
+  splitCandidates.reserve(requestedSize);
+  
+  PolygonContainer::const_iterator it, it_end = polys.end();
+
+  AxisAlignedBox3 box;
+
+  // insert all queries 
+  for(it = polys.begin(); it != it_end; ++ it)
+  {
+	  box.Initialize();
+	  box.Include(*(*it));
+	  
+	  splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
+      splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
+  }
+  
+  stable_sort(splitCandidates.begin(), splitCandidates.end());
+}
+
+
+float BspTree::BestCostRatio(const PolygonContainer &polys,
+							 const AxisAlignedBox3 &box,
+							 const int axis,
+							 float &position,
+							 int &objectsBack,
+							 int &objectsFront) const
+{
+	vector<SortableEntry> splitCandidates;
+
+	SortSplitCandidates(polys, axis, splitCandidates);
+	
+	// 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
+	
+	int objectsLeft = 0, objectsRight = (int)polys.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 = 1e20f;
+	vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
+
+	for(ci = splitCandidates.begin(); ci != ci_end; ++ ci) 
+	{
+		switch ((*ci).type) 
+		{
+			case SortableEntry::POLY_MIN:
+				++ objectsLeft;
+				break;
+			case SortableEntry::POLY_MAX:
+			    -- objectsRight;
+				break;
+			default:
+				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 = objectsLeft * lbox.SurfaceArea() + 
+						objectsRight * rbox.SurfaceArea();
+      
+			if (sum < minSum) 
+			{
+				minSum = sum;
+				position = (*ci).value;
+
+				objectsBack = objectsLeft;
+				objectsFront = objectsRight;
+			}
+		}
+	}
+  
+	float oldCost = (float)polys.size();
+	float newCost = mCt_div_ci + minSum / boxArea;
+	float ratio = newCost / oldCost;
+
+
+#if 0
+  Debug << "====================" << endl;
+  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
+      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
+#endif
+  return ratio;
+}
+
+bool BspTree::SelectAxisAlignedPlane(Plane3 &plane, 
+									 const PolygonContainer &polys) const
+{
+	AxisAlignedBox3 box;
+	box.Initialize();
+	
+	// create bounding box of region
+	Polygon3::IncludeInBox(polys, box);
+        
+	int objectsBack = 0, objectsFront = 0;
+	int axis = 0;
+	float costRatio = MAX_FLOAT;
+	Vector3 position;
+
+	//-- area subdivision
+	for (int i = 0; i < 3; ++ i) 
+	{
+		float p = 0;
+		float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
+		
+		if (r < costRatio)
+		{
+			costRatio = r;
+			axis = i;
+			position = p;
+		}
+	}
+	
+	if (costRatio >= mMaxCostRatio)
+		return false;
+
+	Vector3 norm(0,0,0); norm[axis] = 1.0f;
+	plane = Plane3(norm, position);
+
+	return true;
+}
+
+Plane3 BspTree::SelectPlane(BspLeaf *leaf, BspTraversalData &data)
+{
+	if (data.mPolygons->empty() && data.mRays->empty())
+	{
+		Debug << "Warning: No autopartition polygon candidate available\n";
+	
+		// return axis aligned split
+		AxisAlignedBox3 box;
+		box.Initialize();
+	
+		// create bounding box of region
+		Polygon3::IncludeInBox(*data.mPolygons, box);
+
+		const int axis = box.Size().DrivingAxis();
+		const Vector3 position = (box.Min()[axis] + box.Max()[axis])*0.5f;
+
+		Vector3 norm(0,0,0); norm[axis] = 1.0f;
+		return Plane3(norm, position);
+	}
+	
+	if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
+		((int)data.mPolygons->size() > mTermMaxPolysForAxisAligned) &&
+		((int)data.mRays->size() > mTermMaxRaysForAxisAligned) &&
+		((mTermMaxObjectsForAxisAligned < 0) || 
+		  (Polygon3::ParentObjectsSize(*data.mPolygons) > mTermMaxObjectsForAxisAligned)))
+	{
+		Plane3 plane;
+		if (SelectAxisAlignedPlane(plane, *data.mPolygons))
+			return plane;
+	}
+
+	// simplest strategy: just take next polygon
+	if (mSplitPlaneStrategy & RANDOM_POLYGON)
+	{
+        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();
+	}
+
+	// use heuristics to find appropriate plane
+	return SelectPlaneHeuristics(leaf, data); 
+}
+
+Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf, BspTraversalData &data)
+{
+	float lowestCost = MAX_FLOAT;
+	Plane3 bestPlane;
+	Plane3 plane;
+
+	int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
+	
+	int candidateIdx = limit;
+
+	for (int i = 0; i < limit; ++ i)
+	{
+		candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
+		
+		Polygon3 *poly = (*data.mPolygons)[candidateIdx];
+		// evaluate current candidate
+		const float candidateCost = 
+			SplitPlaneCost(poly->GetSupportingPlane(), data);
+
+		if (candidateCost < lowestCost)
+		{
+			bestPlane = poly->GetSupportingPlane();
+			lowestCost = candidateCost;
+		}
+	}
+	
+	//Debug << "lowest: " << lowestCost << endl;
+
+	//-- choose candidate planes extracted from rays
+	// we currently use two methods
+	// 1) take 3 ray endpoints, where two are minimum and one a maximum
+	//    point or the other way round
+	// 2) take plane normal as plane normal and the midpoint of the ray.
+	//    PROBLEM: does not resemble any point where visibility is likely to change
+	const BoundedRayContainer *rays = data.mRays;
+
+	for (int i = 0; i < mMaxRayCandidates / 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(normal, pt);
+	
+		const float candidateCost = SplitPlaneCost(plane, data);
+
+		if (candidateCost < lowestCost)
+		{
+			bestPlane = plane;
+			
+			lowestCost = candidateCost;
+		}
+	}
+
+	//Debug << "lowest: " << lowestCost << endl;
+	for (int i = 0; i < mMaxRayCandidates / 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;
+		}
+	}	
+
+#ifdef _DEBUG
+	Debug << "plane lowest cost: " << lowestCost << endl;
+#endif
+	return bestPlane;
+}
+
+int BspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
+{
+	const int candidateIdx = Random(currentIdx --);
+
+	// swap candidates to avoid testing same plane 2 times
+	std::swap(polys[currentIdx], polys[candidateIdx]);
+	
+	return currentIdx;
+	//return Random((int)polys.size());
+}
+
+float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
+							  const PolygonContainer &polys) const
+{
+	float val = 0;
+
+	float sumBalancedPolys = 0;
+	float sumSplits = 0;
+	float sumPolyArea = 0;
+	float sumBalancedViewCells = 0;
+	float sumBlockedRays = 0;
+	float totalBlockedRays = 0;
+	//float totalArea = 0;
+	int totalViewCells = 0;
+
+	// need three unique ids for each type of view cell
+	// for balanced view cells criterium
+	ViewCell::NewMail();
+	const int backId = ViewCell::sMailId;
+	ViewCell::NewMail();
+	const int frontId = ViewCell::sMailId;
+	ViewCell::NewMail();
+	const int frontAndBackId = ViewCell::sMailId;
+
+	PolygonContainer::const_iterator it, it_end = polys.end();
+
+	for (it = polys.begin(); it != it_end; ++ it)
+	{
+		const int classification = (*it)->ClassifyPlane(candidatePlane);
+
+		if (mSplitPlaneStrategy & BALANCED_POLYS)
+			sumBalancedPolys += sBalancedPolysTable[classification];
+		
+		if (mSplitPlaneStrategy & LEAST_SPLITS)
+			sumSplits += sLeastPolySplitsTable[classification];
+
+		if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
+		{
+			if (classification == Polygon3::COINCIDENT)
+				sumPolyArea += (*it)->GetArea();
+			//totalArea += area;
+		}
+		
+		if (mSplitPlaneStrategy & BLOCKED_RAYS)
+		{
+			const float blockedRays = (float)(*it)->mPiercingRays.size();
+		
+			if (classification == Polygon3::COINCIDENT)
+				sumBlockedRays += blockedRays;
+			
+			totalBlockedRays += blockedRays;
+		}
+
+		// assign view cells to back or front according to classificaion
+		if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
+		{
+			MeshInstance *viewCell = (*it)->mParent;
+		
+			// assure that we only count a view cell 
+			// once for the front and once for the back side of the plane
+			if (classification == Polygon3::FRONT_SIDE)
+			{
+				if ((viewCell->mMailbox != frontId) && 
+					(viewCell->mMailbox != frontAndBackId))
+				{
+					sumBalancedViewCells += 1.0;
+
+					if (viewCell->mMailbox != backId)
+						viewCell->mMailbox = frontId;
+					else
+						viewCell->mMailbox = frontAndBackId;
+					
+					++ totalViewCells;
+				}
+			}
+			else if (classification == Polygon3::BACK_SIDE)
+			{
+				if ((viewCell->mMailbox != backId) &&
+				    (viewCell->mMailbox != frontAndBackId))
+				{
+					sumBalancedViewCells -= 1.0;
+
+					if (viewCell->mMailbox != frontId)
+						viewCell->mMailbox = backId;
+					else
+						viewCell->mMailbox = frontAndBackId; 
+
+					++ totalViewCells;
+				}
+			}
+		}
+	}
+
+	const float polysSize = (float)polys.size() + Limits::Small;
+
+	// all values should be approx. between 0 and 1 so they can be combined
+	// and scaled with the factors according to their importance
+	if (mSplitPlaneStrategy & BALANCED_POLYS)
+		val += mBalancedPolysFactor * fabs(sumBalancedPolys) / polysSize;
+	
+	if (mSplitPlaneStrategy & LEAST_SPLITS)  
+		val += mLeastSplitsFactor * sumSplits / polysSize;
+
+	if (mSplitPlaneStrategy & LARGEST_POLY_AREA) 
+		// HACK: polys.size should be total area so scaling is between 0 and 1
+		val += mLargestPolyAreaFactor * (float)polys.size() / sumPolyArea;
+
+	if (mSplitPlaneStrategy & BLOCKED_RAYS)
+		if (totalBlockedRays != 0)
+			val += mBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;
+
+	if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
+		val += mBalancedViewCellsFactor * fabs(sumBalancedViewCells) / 
+			((float)totalViewCells + Limits::Small);
+	
+	return val;
+}
+
+bool BspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
+{
+	maxT = 1e6;
+	minT = 0;
+	
+	// test with tree bounding box
+	if (!mBox.GetMinMaxT(ray, &minT, &maxT))
+		return false;
+
+	if (minT < 0) // start ray from origin
+		minT = 0;
+
+	// bound ray or line segment
+	if ((ray.GetType() == Ray::LOCAL_RAY) && 
+	    !ray.intersections.empty() && 
+		(ray.intersections[0].mT <= maxT))
+	{
+		maxT = ray.intersections[0].mT;
+	}
+
+	return true;
+}
+
+float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
+							  const BoundedRayContainer &rays,
+							  const int pvs,
+							  const float area,
+							  const BspNodeGeometry &cell) const
+{
+	float val = 0;
+
+	float sumBalancedRays = 0;
+	float sumRaySplits = 0;
+	
+	int backId = 0;
+	int frontId = 0;
+	int frontAndBackId = 0;
+
+	int frontPvs = 0;
+	int backPvs = 0;
+
+	// probability that view point lies in child
+	float pOverall = 0;
+	float pFront = 0;
+	float pBack = 0;
+
+	if (mSplitPlaneStrategy & PVS)
+	{
+		// create three unique ids for pvs heuristics
+		Intersectable::NewMail(); backId = ViewCell::sMailId;
+		Intersectable::NewMail(); frontId = ViewCell::sMailId;
+		Intersectable::NewMail(); frontAndBackId = ViewCell::sMailId;
+
+		if (mPvsUseArea) // use front and back cell areas to approximate volume
+		{	
+			// construct child geometry with regard to the candidate split plane
+			BspNodeGeometry frontCell;
+			BspNodeGeometry backCell;
+
+			cell.SplitGeometry(frontCell, backCell, *this, candidatePlane);
+		
+			pFront = frontCell.GetArea();
+			pBack = backCell.GetArea();
+
+			pOverall = area;
+		}
+	}
+			
+	BoundedRayContainer::const_iterator rit, rit_end = rays.end();
+
+	for (rit = rays.begin(); rit != rays.end(); ++ rit)
+	{
+		Ray *ray = (*rit)->mRay;
+		const float minT = (*rit)->mMinT;
+		const float maxT = (*rit)->mMaxT;
+
+		Vector3 entP, extP;
+
+		const int cf = 
+			ray->ClassifyPlane(candidatePlane, minT, maxT, entP, extP);
+
+		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+		{
+			sumBalancedRays += sBalancedRaysTable[cf];
+		}
+		
+		if (mSplitPlaneStrategy & BALANCED_RAYS)
+		{
+			sumRaySplits += sLeastRaySplitsTable[cf];
+		}
+
+		if (mSplitPlaneStrategy & PVS)
+		{
+			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, frontPvs, backPvs, 
+					   cf, frontId, backId, frontAndBackId);
+			}
+
+			// the source object in the origin of the ray
+			if (ray->sourceObject.mObject)
+			{
+				IncPvs(*ray->sourceObject.mObject, frontPvs, backPvs, 
+					   cf, frontId, backId, frontAndBackId);
+			}
+
+			if (!mPvsUseArea) // use front and back cell areas to approximate volume
+			{
+				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;
+				}
+			}
+		}
+	}
+
+	const float raysSize = (float)rays.size() + Limits::Small;
+	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+		val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
+
+	if (mSplitPlaneStrategy & BALANCED_RAYS)
+		val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
+
+	float denom = pOverall * (float)pvs * 2.0f + Limits::Small;
+	if ((mSplitPlaneStrategy & PVS) && area && pvs)
+	{
+		val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
+
+		// give penalty to unbalanced split
+		if (0)
+		if (((pFront * 0.2 + Limits::Small) > pBack) || (pFront < (pBack * 0.2 + Limits::Small)))
+			val += 0.5;
+	}
+
+#ifdef _DEBUG
+	Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
+		  << " frontpvs: " << frontPvs << " pFront: " << pFront 
+		  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
+#endif
+	return val;
+}
+
+void BspTree::IncPvs(Intersectable &obj,
+					 int &frontPvs,
+					 int &backPvs,
+					 const int cf, 
+					 const int frontId, 
+					 const int backId, 
+					 const int frontAndBackId) const
+{
+	// TODO: does this really belong to no pvs?
+	//if (cf == Ray::COINCIDENT) return;
+
+	if (cf == Ray::FRONT)
+	{
+		if ((obj.mMailbox != frontId) && 
+			(obj.mMailbox != frontAndBackId))
+		{
+			++ frontPvs;
+
+			if (obj.mMailbox != backId)
+				obj.mMailbox = frontId;
+			else
+				obj.mMailbox = frontAndBackId;					
+		}
+	}
+	else if (cf == Ray::BACK)
+	{
+		if ((obj.mMailbox != backId) &&
+			(obj.mMailbox != frontAndBackId))
+		{
+			++ backPvs;
+
+			if (obj.mMailbox != frontId)
+				obj.mMailbox = backId;
+			else
+				obj.mMailbox = frontAndBackId; 
+		}
+	}
+	// object belongs to both PVS
+	else if ((cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT) ||(cf == Ray::COINCIDENT))
+	{
+		if (obj.mMailbox !=  frontAndBackId)
+		{
+			if (obj.mMailbox != frontId)
+				++ frontPvs;
+			if (obj.mMailbox != backId)
+				++ backPvs;
+		
+			obj.mMailbox = frontAndBackId;
+		}
+	}
+}
+
+float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
+							  BspTraversalData &data) const
+{
+	float val = 0;
+
+	if (mSplitPlaneStrategy & VERTICAL_AXIS)
+	{
+		Vector3 tinyAxis(0,0,0); tinyAxis[mBox.Size().TinyAxis()] = 1.0f;
+		// we put a penalty on the dot product between the "tiny" vertical axis
+		// and the split plane axis
+		val += mVerticalSplitsFactor * 
+			   fabs(DotProd(candidatePlane.mNormal, tinyAxis));
+	}
+
+	// the following criteria loop over all polygons to find the cost value
+	if ((mSplitPlaneStrategy & BALANCED_POLYS)      ||
+		(mSplitPlaneStrategy & LEAST_SPLITS)        ||
+		(mSplitPlaneStrategy & LARGEST_POLY_AREA)   ||
+		(mSplitPlaneStrategy & BALANCED_VIEW_CELLS) ||
+		(mSplitPlaneStrategy & BLOCKED_RAYS))
+	{
+		val += SplitPlaneCost(candidatePlane, *data.mPolygons);
+	}
+
+	// the following criteria loop over all rays to find the cost value
+	if ((mSplitPlaneStrategy & BALANCED_RAYS)      ||
+		(mSplitPlaneStrategy & LEAST_RAY_SPLITS)   ||
+		(mSplitPlaneStrategy & PVS))
+	{
+		val += SplitPlaneCost(candidatePlane, *data.mRays, data.mPvs, 
+							  data.mArea, *data.mGeometry);
+	}
+
+	// return linear combination of the sums
+	return val;
 }
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp	(revision 421)
@@ -23,7 +23,6 @@
   environment->Parse(argc, argv, USE_EXE_PATH);
   MeshKdTree::ParseEnvironment();
-  
-  VspKdTree::ParseEnvironment();
-
+  BspTree::ParseEnvironment();
+ 
   char buff[128];
   
@@ -63,7 +62,7 @@
   // a certain number of rays collected
   if (ViewCell::sHierarchy == ViewCell::BSP &&
-	  !(BspTree::mConstructionMethod == BspTree::FROM_RAYS)) 
+	  !(BspTree::sConstructionMethod == BspTree::FROM_RAYS)) 
   {
-	  if (BspTree::mConstructionMethod == BspTree::FROM_INPUT_VIEW_CELLS)
+	  if (BspTree::sConstructionMethod == BspTree::FROM_INPUT_VIEW_CELLS)
 	  {
 		  // view cells input file name
Index: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h	(revision 421)
@@ -421,29 +421,4 @@
 	float GetAvgPvsSize();
 
-	/** Parses the environment and stores the global VspKd tree parameters
-	*/
-	static void ParseEnvironment();
-
-	/////////////////////////////
-	// Construction parameters
-
-	/// max depth of the tree
-	static int sTermMaxDepth;
-
-	/// minimal ratio of the volume of the cell and the query volume
-	static float sTermMinSize;
-
-	/// minimal pvs per node to still get subdivided
-	static int sTermMinPvs;
-
-	/// minimal ray number per node to still get subdivided
-	static int sTermMinRays;
-	
-	/// maximal cost ration to subdivide a node
-	static float sTermMaxCostRatio;
-	
-	/// maximal contribution per ray to subdivide the node
-	static float sTermMaxRayContribution;
-
 	/** Returns memory usage in MB.
 	*/
@@ -533,4 +508,7 @@
 	int ReleaseMemory(const int time);
 
+	bool TerminationCriteriaMet(const VspKdTreeLeaf *leaf, 
+								const AxisAlignedBox3 &box) const;
+
 protected:
 	/////////////////////////////
@@ -575,5 +553,28 @@
 	// reusable array of split candidates
 	vector<SortableEntry> *splitCandidates;
-	
+
+	
+	/////////////////////////////
+	// Construction parameters
+
+	/// max depth of the tree
+	int mTermMaxDepth;
+
+	/// minimal ratio of the volume of the cell and the query volume
+	float mTermMinSize;
+
+	/// minimal pvs per node to still get subdivided
+	int mTermMinPvs;
+
+	/// minimal ray number per node to still get subdivided
+	int mTermMinRays;
+	
+	/// maximal cost ration to subdivide a node
+	float mTermMaxCostRatio;
+	
+	/// maximal contribution per ray to subdivide the node
+	float mTermMaxRayContribution;
+
+
 	/////////////////////////////
 	VspKdStatistics mStat;	
Index: trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 420)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 421)
@@ -61,5 +61,5 @@
 	ObjectContainer objects;
 	
-	switch (mBspTree->mConstructionMethod)
+	switch (BspTree::sConstructionMethod)
 	{
 	case BspTree::FROM_INPUT_VIEW_CELLS:
@@ -779,5 +779,5 @@
 	if (!mBspTree)
 	{
-		if ((BspTree::mConstructionMethod == BspTree::FROM_RAYS) &&
+		if ((BspTree::sConstructionMethod == BspTree::FROM_RAYS) &&
 			((int)mSampleRays.size() < mBspConstructionSamples))
 		{
@@ -889,5 +889,5 @@
 					outRays.push_back(mSampleRays[i]);
 			}
-			if (BspTree::mConstructionMethod == BspTree::FROM_RAYS)
+			if (BspTree::sConstructionMethod == BspTree::FROM_RAYS)
 			{
 				// export rays 
