Index: /trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj	(revision 439)
+++ /trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj	(revision 448)
@@ -255,4 +255,16 @@
 			</File>
 			<File
+				RelativePath="..\src\RssPreprocessor.cpp">
+			</File>
+			<File
+				RelativePath="..\src\RssPreprocessor.h">
+			</File>
+			<File
+				RelativePath="..\src\RssTree.cpp">
+			</File>
+			<File
+				RelativePath="..\src\RssTree.h">
+			</File>
+			<File
 				RelativePath="..\src\SamplingPreprocessor.cpp">
 			</File>
Index: /trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 445)
+++ /trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 448)
@@ -19,6 +19,6 @@
 
 Preprocessor {
-	type sampling
-#	type vss
+#	type sampling
+	type vss
 }
 
@@ -94,5 +94,5 @@
 
 Sampling {
-	totalSamples 600000
+	totalSamples 300000
 	samplesPerPass  3
 }
@@ -101,7 +101,7 @@
 	loadFromFile false
 	#type kdTree
-	#type vspTree
+	type vspKdTree
 	#type bspTree
-	type vspBspTree
+	#type vspBspTree
 	
 	#type sceneDependent
@@ -131,6 +131,6 @@
 BspTree {
 	Construction {
-		samples 200000
-		sideTolerance 0.005
+		samples 50000
+		epsilon 0.005
 	}
 
@@ -193,5 +193,5 @@
 		minRays 200
 		minPolygons -1
-		maxDepth 50
+		maxDepth 40
 		minPvs 100
 		minArea 0.01
@@ -256,5 +256,5 @@
 
 	Construction {
-		samples 500000
+		samples 100000
 	}
 	
@@ -278,6 +278,6 @@
 VspBspTree {
 	Construction {
-		samples 200000
-		sideTolerance 0.005
+		samples 100000
+		epsilon 0.005
 	}
 
@@ -309,5 +309,5 @@
 		minRays 200
 		minPolygons -1
-		maxDepth 50
+		maxDepth 40
 		minPvs 100
 		minArea 0.01
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.h	(revision 445)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.h	(revision 448)
@@ -9,4 +9,5 @@
 #include "VssRay.h"
 #include "RayInfo.h"
+#include "ViewCellBsp.h"
 
 class ViewCell;
@@ -14,318 +15,16 @@
 class Plane3;
 class VspBspTree;  
-class VspBspInterior;
-class VspBspNode;
+class BspInterior;
+class BspNode;
 class AxisAlignedBox3;
 class Ray;
 
-class VspBspNodeGeometry
-{
-public:
-	VspBspNodeGeometry()
-	{};  
-
-	~VspBspNodeGeometry();
-
-	float GetArea() const;
-	
-	/** Computes new cell based on the old cell definition and a new split plane
-		@param side indicates which side of the halfspace 
-	*/
-	void SplitGeometry(VspBspNodeGeometry &front, 
-					   VspBspNodeGeometry &back,
-					   const VspBspTree &tree,
-					   const Plane3 &splitPlane) const;
-
-	Polygon3 *SplitPolygon(Polygon3 *poly, const VspBspTree &tree) const;
-
-	PolygonContainer mPolys;
-};
-
-/** Data structure used for optimized ray casting.
+/*class BspNodeGeometry;
+class BspTreeStatistics;
+class BspViewCellsStatistics;
+class BspNode;
+class BspLeaf;
+class BspInterior;
 */
-struct VspBspRayTraversalData 
-{
-    VspBspNode *mNode;
-    Vector3 mExitPoint;
-    float mMaxT;
-    
-    VspBspRayTraversalData() {}
-
-    VspBspRayTraversalData(VspBspNode *n, const Vector3 &extp, const float maxt):
-    mNode(n), mExitPoint(extp), mMaxT(maxt) 
-	{}
-};
-
-class VspBspTreeStatistics: public StatisticsBase
-{
-public:
-	// total number of nodes
-	int nodes;
-	// number of splits
-	int splits;
-	// totals number of rays
-	int rays;
-	// maximal reached depth
-	int maxDepth;
-	// minimal depth
-	int minDepth;
-	
-	// max depth nodes
-	int maxDepthNodes;
-	// minimum depth nodes
-	int minDepthNodes;
-	// max depth nodes
-	int minPvsNodes;
-	// nodes with minimum PVS
-	int minRaysNodes;
-	// max ray contribution nodes
-	int maxRayContribNodes;
-	// minimum area nodes
-	int minAreaNodes;
-
-	// max number of rays per node
-	int maxObjectRefs;
-	// accumulated depth (used to compute average)
-	int accumDepth;
-	// number of initial polygons
-	int polys;
-	/// samples contributing to pvs
-	int contributingSamples;
-	/// sample contributions to pvs
-	int sampleContributions;
-	/// largest pvs
-	int largestPvs;
-
-	// Constructor
-	VspBspTreeStatistics() 
-	{
-		Reset();
-	}
-
-	int Nodes() const {return nodes;}
-	int Interior() const { return nodes / 2; }
-	int Leaves() const { return (nodes / 2) + 1; }
-	
-	// TODO: computation wrong
-	double AvgDepth() const { return accumDepth / (double)Leaves();}; 
-  
-	void Reset() 
-	{
-		nodes = 0;
-		splits = 0;
-		
-		maxDepth = 0;
-		minDepth = 99999;
-		polys = 0;
-		accumDepth = 0;
-        
-		maxDepthNodes = 0;
-		minPvsNodes = 0;
-		minRaysNodes = 0;
-		maxRayContribNodes = 0;
-		minAreaNodes = 0;
-
-		contributingSamples = 0;
-		sampleContributions = 0;
-	}
-
-	void Print(ostream &app) const;
-
-	friend ostream &operator<<(ostream &s, const VspBspTreeStatistics &stat) 
-	{
-		stat.Print(s);
-		return s;
-	} 
-};
-
-class VspBspViewCellsStatistics: public StatisticsBase
-{
-public:
-
-	/// number of view cells
-	int viewCells;
-
-	/// size of the PVS
-	int pvs;
-  
-	/// largest PVS of all view cells
-	int maxPvs;
-
-	/// smallest PVS of all view cells
-	int minPvs;
-
-	/// view cells with empty PVS
-	int emptyPvs;
-
-	/// number of bsp leaves covering the view space
-	int bspLeaves;
-
-	/// largest number of leaves covered by one view cell  
-	int maxVspBspLeaves;
-
-    // Constructor
-	VspBspViewCellsStatistics() 
-	{
-		Reset();
-	}
-
-	double AvgVspBspLeaves() const {return (double)bspLeaves / (double)viewCells;};
-	double AvgPvs() const {return (double)pvs / (double)viewCells;};
-  
-	void Reset() 
-	{
-		viewCells = 0;	
-		pvs = 0;  
-		maxPvs = 0;
-
-		minPvs = 999999;
-		emptyPvs = 0;
-		bspLeaves = 0;
-		maxVspBspLeaves = 0;
-	}
-
-	void Print(ostream &app) const;
-
-	friend ostream &operator<<(ostream &s, const VspBspViewCellsStatistics &stat) 
-	{
-		stat.Print(s);
-		return s;
-	}  
-};
-
-/**
-    VspBspNode abstract class serving for interior and leaf node implementation
-*/
-class VspBspNode 
-{
-	friend class VspBspTree;
-
-public:
-	VspBspNode();
-	virtual ~VspBspNode(){};
-	VspBspNode(VspBspInterior *parent);
-
-	/** Determines whether this node is a leaf or not
-	@return true if leaf
-	*/
-	virtual bool IsLeaf() const = 0;
-
-	/** Determines whether this node is a root
-	@return true if root
-	*/
-	virtual bool IsRoot() const;
-
-	/** Returns parent node.
-	*/
-	VspBspInterior *GetParent();
-
-	/** Sets parent node.
-	*/
-	void SetParent(VspBspInterior *parent);
-
-	
-	static int sMailId;
-	int mMailbox;
-  
-	void Mail() { mMailbox = sMailId; }
-	static void NewMail() { ++ sMailId; }
-	bool Mailed() const { return mMailbox == sMailId; }
-
-protected:
-
-	/// parent of this node
-	VspBspInterior *mParent;
-};
-
-/** BSP interior node implementation 
-*/
-class VspBspInterior : public VspBspNode 
-{
-	friend class VspBspTree;
-public:
-	/** Standard contructor taking split plane as argument.
-	*/
-	VspBspInterior(const Plane3 &plane);
-	~VspBspInterior();
-	/** @return false since it is an interior node 
-	*/
-	bool IsLeaf() const;
-
-	VspBspNode *GetBack();
-	VspBspNode *GetFront();
-
-	Plane3 *GetPlane();
-
-	void ReplaceChildLink(VspBspNode *oldChild, VspBspNode *newChild);
-	void SetupChildLinks(VspBspNode *b, VspBspNode *f);
-
-	/** Splits polygons with respect to the split plane.
-		@param polys the polygons to be split. the polygons are consumed and
-			   distributed to the containers frontPolys, backPolys, coincident.
-		@param frontPolys returns the polygons in the front of the split plane
-		@param backPolys returns the polygons in the back of the split plane
-		@param coincident returns the polygons coincident to the split plane
-
-		@returns the number of splits	
-	*/
-	int SplitPolygons(PolygonContainer &polys, 
-					  PolygonContainer &frontPolys, 
-					  PolygonContainer &backPolys, 
-					  PolygonContainer &coincident);
-
-	friend ostream &operator<<(ostream &s, const VspBspInterior &A)
-	{
-		return s << A.mPlane;
-	}
-
-protected:
-
-	/// Splitting plane corresponding to this node
-	Plane3 mPlane;
-	/// back node
-	VspBspNode *mBack;
-	/// front node
-	VspBspNode *mFront;
-};
-
-/** BSP leaf node implementation.
-*/
-class VspBspLeaf : public VspBspNode 
-{
-	friend class VspBspTree;
-
-public:
-	VspBspLeaf();
-	VspBspLeaf(BspViewCell *viewCell);
-	VspBspLeaf(VspBspInterior *parent);
-	VspBspLeaf(VspBspInterior *parent, BspViewCell *viewCell);
-
-	/** @return true since it is an interior node 
-	*/
-	bool IsLeaf() const;
-	
-	/** Returns pointer of view cell.
-	*/
-	BspViewCell *GetViewCell() const;
-
-	/** Sets pointer to view cell.
-	*/
-	void SetViewCell(BspViewCell *viewCell);
-
-	/** Adds ray sample contributions to the PVS.
-		@param sampleContributions the number contributions of the samples
-		@param contributingSampels the number of contributing rays
-		
-	*/
-	void AddToPvs(const RayInfoContainer &rays, 
-				  int &sampleContributions,
-				  int &contributingSamples, 
-				  bool storeLeavesWithRays = false);
-
-protected:
-
-	/// if NULL this does not correspond to feasible viewcell
-	BspViewCell *mViewCell;
-};
 
 /**
@@ -345,18 +44,17 @@
 	{  
 		/// the current node
-		VspBspNode *mNode;
+		BspNode *mNode;
 		/// polygonal data for splitting
 		PolygonContainer *mPolygons;
 		/// current depth
 		int mDepth;
-		/// the view cell associated with this subdivsion
-		ViewCell *mViewCell;
+		
 		/// rays piercing this node
 		RayInfoContainer *mRays;
 		/// area of current node
 		float mArea;
-		/// geometry of current node induced by split planes
-		VspBspNodeGeometry *mGeometry;
-
+		/// geometry of node as induced by planes
+		BspNodeGeometry *mGeometry;
+		
 		/// pvs size
 		int mPvs;
@@ -374,5 +72,4 @@
 		mPolygons(NULL),
 		mDepth(0),
-		mViewCell(NULL),
 		mRays(NULL),
 		mPvs(0),
@@ -381,20 +78,31 @@
 		{}
 		
-		VspBspTraversalData(VspBspNode *node, 
+		VspBspTraversalData(BspNode *node, 
 						 PolygonContainer *polys, 
 						 const int depth, 
-						 ViewCell *viewCell,
 						 RayInfoContainer *rays,
 						 int pvs,
 						 float area,
-						 VspBspNodeGeometry *cell): 
+						 BspNodeGeometry *geom): 
 		mNode(node), 
 		mPolygons(polys), 
 		mDepth(depth), 
-		mViewCell(viewCell),
 		mRays(rays),
 		mPvs(pvs),
 		mArea(area),
-		mGeometry(cell)
+		mGeometry(geom)
+		{}
+
+		VspBspTraversalData(PolygonContainer *polys, 
+							const int depth, 
+							RayInfoContainer *rays,
+							BspNodeGeometry *geom): 
+		mNode(NULL), 
+		mPolygons(polys), 
+		mDepth(depth), 
+		mRays(rays),
+		mPvs(0),
+		mArea(0),
+		mGeometry(geom)
 		{}
     };
@@ -406,7 +114,11 @@
 	VspBspTree();
 
+	/** Default destructor.
+	*/
 	~VspBspTree();
 
-	const VspBspTreeStatistics &GetStatistics() const; 
+	/** Returns BSP Tree statistics.
+	*/
+	const BspTreeStatistics &GetStatistics() const; 
   
 
@@ -420,5 +132,5 @@
 	/** Returns list of BSP leaves.
 	*/
-	void CollectLeaves(vector<VspBspLeaf *> &leaves) const;
+	void CollectLeaves(vector<BspLeaf *> &leaves) const;
 
 	/** Returns box which bounds the whole tree.
@@ -428,5 +140,5 @@
 	/** Returns root of BSP tree.
 	*/
-	VspBspNode *GetRoot() const;
+	BspNode *GetRoot() const;
 
 	/** Exports VspBsp tree to file.
@@ -450,9 +162,10 @@
 	/** Returns statistics.
 	*/
-	VspBspTreeStatistics &GetStat();
+	BspTreeStatistics &GetStat();
 
 	/** finds neighbouring leaves of this tree node.
 	*/
-	int FindNeighbors(VspBspNode *n, vector<VspBspLeaf *> &neighbors, 
+	int FindNeighbors(BspNode *n, 
+					  vector<BspLeaf *> &neighbors, 
 					  const bool onlyUnmailed) const;
 
@@ -460,5 +173,5 @@
 		leading to this node.
 	*/
-	void ConstructGeometry(VspBspNode *n, PolygonContainer &cell) const;
+	void ConstructGeometry(BspNode *n, PolygonContainer &cell) const;
 
 	/** Constructs geometry associated with the half space intersections 
@@ -469,19 +182,19 @@
 	/** Construct geometry and stores it in a geometry node container.
 	*/
-	void ConstructGeometry(VspBspNode *n, VspBspNodeGeometry &cell) const;
+	void ConstructGeometry(BspNode *n, BspNodeGeometry &cell) const;
 
 	/** Returns random leaf of BSP tree.
 		@param halfspace defines the halfspace from which the leaf is taken.
 	*/
-	VspBspLeaf *GetRandomLeaf(const Plane3 &halfspace);
+	BspLeaf *GetRandomLeaf(const Plane3 &halfspace);
 
 	/** Returns random leaf of BSP tree.
 		@param onlyUnmailed if only unmailed leaves should be returned.
 	*/
-	VspBspLeaf *GetRandomLeaf(const bool onlyUnmailed = false);
+	BspLeaf *GetRandomLeaf(const bool onlyUnmailed = false);
 
 	/** Traverses tree and counts all view cells as well as their PVS size.
 	*/
-	void EvaluateViewCellsStats(VspBspViewCellsStatistics &stat) const;
+	void EvaluateViewCellsStats(BspViewCellsStatistics &stat) const;
 
 
@@ -489,4 +202,8 @@
 	*/
 	BspViewCell *GetRootCell() const;
+
+	/** Returns epsilon of this tree.
+	*/
+	float GetEpsilon() const;
 
 protected:
@@ -521,11 +238,12 @@
 		@returns new root of the subtree
 	*/
-	VspBspNode *Subdivide(VspBspTraversalStack &tStack, VspBspTraversalData &tData);
-
-	/** Constructs the tree from the given list of polygons and rays.
+	BspNode *Subdivide(VspBspTraversalStack &tStack, 
+					   VspBspTraversalData &tData);
+
+	/** Constructs the tree from the given traversal data.
 		@param polys stores set of polygons on which subdivision may be based
 		@param rays storesset of rays on which subdivision may be based
 	*/
-	void Construct(PolygonContainer *polys, RayInfoContainer *rays);
+	void Construct(const PolygonContainer &polys, RayInfoContainer *rays);
 
 	/** Selects the best possible splitting plane. 
@@ -536,5 +254,5 @@
 		@returns the split plane
 	*/
-	Plane3 SelectPlane(VspBspLeaf *leaf, 
+	Plane3 SelectPlane(BspLeaf *leaf, 
 					   VspBspTraversalData &data);
 
@@ -564,5 +282,5 @@
 	*/
 
-	VspBspInterior *SubdivideNode(VspBspTraversalData &tData,
+	BspInterior *SubdivideNode(VspBspTraversalData &tData,
 							   VspBspTraversalData &frontData,
 							   VspBspTraversalData &backData,
@@ -575,5 +293,5 @@
 		@param rays bundle of rays on which the split can be based
 	*/
-	Plane3 SelectPlaneHeuristics(VspBspLeaf *leaf,
+	Plane3 SelectPlaneHeuristics(BspLeaf *leaf,
 								 VspBspTraversalData &data);
 
@@ -639,9 +357,4 @@
 	bool SelectAxisAlignedPlane(Plane3 &plane, const PolygonContainer &polys) const;
 
-	/** Bounds ray and returns minT and maxT.
-		@returns true if ray hits BSP tree bounding box
-	*/
-	bool BoundRay(const Ray &ray, float &minT, float &maxT) const;
-
 	/** Subdivides the rays into front and back rays according to the split plane.
 		
@@ -662,5 +375,5 @@
 	/** Extracts the split planes representing the space bounded by node n.
 	*/
-	void ExtractHalfSpaces(VspBspNode *n, vector<Plane3> &halfSpaces) const;
+	void ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const;
 
 	/** Adds the object to the pvs of the front and back leaf with a given classification.
@@ -686,10 +399,35 @@
 	float AccumulatedRayLength(const RayInfoContainer &rays) const;
 
-
+	/** Splits polygons with respect to the split plane.
+
+		@param plane the split plane
+		@param polys the polygons to be split. the polygons are consumed and
+			   distributed to the containers frontPolys, backPolys, coincident.
+		@param frontPolys returns the polygons in the front of the split plane
+		@param backPolys returns the polygons in the back of the split plane
+		@param coincident returns the polygons coincident to the split plane
+
+		@returns the number of splits	
+	*/
+	int SplitPolygons(const Plane3 &plane,
+					  PolygonContainer &polys, 
+					  PolygonContainer &frontPolys, 
+					  PolygonContainer &backPolys, 
+					  PolygonContainer &coincident) const;
+
+	/** Adds ray sample contributions to the PVS.
+		@param sampleContributions the number contributions of the samples
+		@param contributingSampels the number of contributing rays
+		
+	*/
+	void AddToPvs(BspLeaf *leaf,
+				  const RayInfoContainer &rays, 
+				  int &sampleContributions,
+				  int &contributingSamples);
 
 	/// Pointer to the root of the tree
-	VspBspNode *mRoot;
-		
-	VspBspTreeStatistics mStat;
+	BspNode *mRoot;
+		
+	BspTreeStatistics mStat;
 
 	/// Strategies for choosing next split plane.
@@ -749,4 +487,8 @@
 	bool mPvsUseArea;
 
+	float mEpsilon;
+
+	int mMaxTests;
+
 private:
 	
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp	(revision 440)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp	(revision 448)
@@ -459,19 +459,36 @@
 	
 	mStat.nodes = 1;
-  	
-	mBox.Initialize();
-	
-	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
-	{
-		leaf->AddRay(RayInfo(*ri));
-
-		mBox.Include((*ri)->GetOrigin());
-		mBox.Include((*ri)->GetTermination()); 
-	}
-	
+  	mBox.Initialize();
+
+	//-- compute bounding box
 	if (forcedBoundingBox)
 		mBox = *forcedBoundingBox;
+	else
+		for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
+		{
+			mBox.Include((*ri)->GetOrigin());
+			mBox.Include((*ri)->GetTermination()); 
+		}
 
 	Debug << "bbox = " << mBox << endl;
+
+	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
+	{
+		//leaf->AddRay(RayInfo(*ri));
+		VssRay *ray = *ri;
+		
+		float minT, maxT;
+
+		// TODO: not very efficient to implictly cast between rays types ...
+		if (mBox.GetRaySegment(*ray, minT, maxT))
+		{
+			float len = ray->Length();
+			
+			if (!len) 
+				len = Limits::Small;
+		
+			leaf->AddRay(RayInfo(ray, minT / len, maxT / len)); 
+		}
+	}
 	
 	mStat.rays = (int)leaf->mRays.size();
@@ -1317,5 +1334,5 @@
 									 stack<RayTraversalData> &tstack)
 {
-	VspKdTreeInterior *in =  (VspKdTreeInterior *) data.mNode;
+	VspKdTreeInterior *in = dynamic_cast<VspKdTreeInterior *>(data.mNode);
   
 	if (in->mAxis <= VspKdTreeNode::SPLIT_Z) 
@@ -1688,5 +1705,5 @@
     
 	AxisAlignedBox3 box(node->mParent->mBox);
-	if (node->mParent->mFront == node)
+	if (node->mParent->GetFront() == node)
 		box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
 	else
@@ -1763,7 +1780,7 @@
 			VspKdTreeInterior *interior = dynamic_cast<VspKdTreeInterior *>(node);
 
-			nodeStack.push(interior->mBack);
-			nodeStack.push(interior->mFront);
-		}
-	}
-}
+			nodeStack.push(interior->GetBack());
+			nodeStack.push(interior->GetFront());
+		}
+	}
+}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h	(revision 436)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h	(revision 448)
@@ -50,8 +50,10 @@
 		@param front returns the front the front polygon
 		@param back returns the back polygon
+		@param epsilon epsilon where two points are considered equal
 	*/
 	void Split(const Plane3 &partition, 
 			   Polygon3 &front, 
-			   Polygon3 &back);
+			   Polygon3 &back,
+			   const float epsilon);// = Limits::Small);
 
 	/** Returns the area of this polygon.
@@ -62,16 +64,17 @@
 	    @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
 	*/
-	int ClassifyPlane(const Plane3 &plane) const;
-
+	//int ClassifyPlane(const Plane3 &plane) const;
 	
 	/** Classify polygon with respect to the plane.
-	    @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
+	    @param epsilon tolerance value
+		@returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
+
 	*/
-	int ClassifyPlane(const Plane3 &plane, const float tolerance) const;
+	int ClassifyPlane(const Plane3 &plane, const float epsilon) const;
 
 	/** Side of the polygon with respect to the plane.
 		@returns 1 if on front side, -1 if on back side, 0 else.
 	*/
-	int Side(const Plane3 &plane) const;
+	int Side(const Plane3 &plane, const float epsilon) const;
 
 	/** Scales the polygon about its center
@@ -85,5 +88,5 @@
 		@returns true if polygon is valid.
 	*/
-	bool Valid() const; 
+	bool Valid(const float epsilon) const; 
 
 	/** Returns the surface normal.
@@ -127,5 +130,7 @@
 	    @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
 	*/
-	static int ClassifyPlane(const PolygonContainer &polys, const Plane3 &plane);
+	static int ClassifyPlane(const PolygonContainer &polys, 
+							 const Plane3 &plane,
+							 const float epsilon);
 
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 444)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 448)
@@ -7,4 +7,6 @@
 #include <stack>
 #include "Statistics.h"
+#include "VssRay.h"
+
 
 class ViewCell;
@@ -27,13 +29,14 @@
 	float GetArea() const;
 	
-	/** Computes new cell based on the old cell definition and a new split plane
-		@param side indicates which side of the halfspace 
-	*/
-	void SplitGeometry(BspNodeGeometry &front, 
+	/** Computes new front and back geometry based on the old cell 
+		geometry and a new split plane
+	*/
+	void SplitGeometry(BspNodeGeometry &front,
 					   BspNodeGeometry &back,
-					   const BspTree &tree,
-					   const Plane3 &splitPlane) const;
-
-	Polygon3 *SplitPolygon(Polygon3 *poly, const BspTree &tree) const;
+					   const Plane3 &splitPlane,
+                       const AxisAlignedBox3 &box,
+					   const float epsilon) const;
+
+	Polygon3 *SplitPolygon(Polygon3 *poly, const float epsilon) const;
 
 	PolygonContainer mPolys;
@@ -270,22 +273,14 @@
 	BspNode *GetFront();
 
-	Plane3 *GetPlane();
-
+	/** Returns split plane.
+	*/
+	Plane3 GetPlane() const;
+
+	/** Replace front or back child with new child.
+	*/
 	void ReplaceChildLink(BspNode *oldChild, BspNode *newChild);
+	/** Replace front and back child.
+	*/
 	void SetupChildLinks(BspNode *b, BspNode *f);
-
-	/** Splits polygons with respect to the split plane.
-		@param polys the polygons to be split. the polygons are consumed and
-			   distributed to the containers frontPolys, backPolys, coincident.
-		@param frontPolys returns the polygons in the front of the split plane
-		@param backPolys returns the polygons in the back of the split plane
-		@param coincident returns the polygons coincident to the split plane
-
-		@returns the number of splits	
-	*/
-	int SplitPolygons(PolygonContainer &polys, 
-					  PolygonContainer &frontPolys, 
-					  PolygonContainer &backPolys, 
-					  PolygonContainer &coincident);
 
 	friend ostream &operator<<(ostream &s, const BspInterior &A)
@@ -298,4 +293,5 @@
 	/// Splitting plane corresponding to this node
 	Plane3 mPlane;
+
 	/// back node
 	BspNode *mBack;
@@ -328,12 +324,5 @@
 	void SetViewCell(BspViewCell *viewCell);
 
-	/** Adds ray sample contributions to the PVS.
-		@param sampleContributions the number contributions of the samples
-		@param contributingSampels the number of contributing rays
-		
-	*/
-	void AddToPvs(const BoundedRayContainer &rays, 
-				  int &sampleContributions,
-				  int &contributingSamples);
+	VssRayContainer mVssRays;
 
 protected:
@@ -516,8 +505,7 @@
 	BspViewCell *GetRootCell() const;
 
-	/** Parses the environment and stores the global BSP tree parameters
-	*/
-	static void ParseEnvironment();
-
+	/** Returns epsilon of this tree.
+	*/
+	float GetEpsilon() const;
 
 protected:
@@ -722,9 +710,4 @@
 	bool SelectAxisAlignedPlane(Plane3 &plane, const PolygonContainer &polys) const;
 
-	/** Bounds ray and returns minT and maxT.
-		@returns true if ray hits BSP tree bounding box
-	*/
-	bool BoundRay(const Ray &ray, float &minT, float &maxT) const;
-
 	/** Subdivides the rays into front and back rays according to the split plane.
 		
@@ -768,4 +751,29 @@
 	*/
 	float AccumulatedRayLength(BoundedRayContainer &rays) const;
+
+	/** Splits polygons with respect to the split plane.
+		@param polys the polygons to be split. the polygons are consumed and
+			   distributed to the containers frontPolys, backPolys, coincident.
+		@param frontPolys returns the polygons in the front of the split plane
+		@param backPolys returns the polygons in the back of the split plane
+		@param coincident returns the polygons coincident to the split plane
+
+		@returns the number of splits	
+	*/
+	int SplitPolygons(const Plane3 &plane,
+					  PolygonContainer &polys, 
+					  PolygonContainer &frontPolys, 
+					  PolygonContainer &backPolys, 
+					  PolygonContainer &coincident) const;
+
+	/** Adds ray sample contributions to the PVS.
+		@param sampleContributions the number contributions of the samples
+		@param contributingSampels the number of contributing rays
+		
+	*/
+	void AddToPvs(BspLeaf *leaf,
+				  const BoundedRayContainer &rays, 
+				  int &sampleContributions,
+				  int &contributingSamples);
 
 	/// Pointer to the root of the tree.
@@ -852,4 +860,7 @@
 	bool mPvsUseArea;
 
+	/// epsilon where two points are still considered equal
+	float mEpsilon;
+
 private:
 	
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h	(revision 441)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h	(revision 448)
@@ -9,4 +9,5 @@
 class ViewCell;
 class KdLeaf;
+class VspBspTree;
 
 class SimulationStatistics: public StatisticsBase
@@ -105,5 +106,29 @@
 	Real RenderPvs(ViewCell &viewCell, const float objRenderTime) const;
 
+private:
 	BspTree *mBspTree;
+};
+
+class VspBspRenderSimulator: public RenderSimulator
+{
+public:
+	VspBspRenderSimulator(VspBspTree *vspBspTree);
+
+	VspBspRenderSimulator(float objRenderCost, 
+						  float vcOverhead, 
+						  float moveSpeed, 
+						  VspBspTree *vspBspTree);
+
+	SimulationStatistics SimulateRendering();
+
+protected:
+	/** Simulates rendering of the pvs of one view cell, with given rendering time for an object.
+		@param viewCell the view cell holding the Pvs
+		@param objRenderTime estimated render time for one object of the Pvs
+	*/
+	Real RenderPvs(ViewCell &viewCell, const float objRenderTime) const;
+
+private:
+	VspBspTree *mVspBspTree;
 };
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h	(revision 387)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h	(revision 448)
@@ -312,5 +312,12 @@
 	  @returns the cross section
   */
-  Polygon3 *CrossSection(const Plane3 &plane);
+  Polygon3 *CrossSection(const Plane3 &plane) const;
+
+  /** Computes minimal and maximal t of ray, including the object intersections.
+	  @returns true if ray hits the bounding box.
+  */
+  bool GetRaySegment(const Ray &ray, 
+				float &minT, 
+				float &maxT) const;
 
 #define __EXTENT_HACK
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h	(revision 448)
@@ -36,8 +36,4 @@
   Vector3(float X) { x = y = z = X; }
   Vector3(const Vector3 &v) { x = v.x; y = v.y; z = v.z; }
-
-  /// the distance where two points are still considered equal
-  static float sDistTolerance;
-  static float sDistToleranceSqrt;
 
   // Functions to get at the vector components
@@ -480,5 +476,5 @@
 EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2)
 {
-  return EpsilonEqualV3(v1,v2,Limits::Small);
+  return EpsilonEqualV3(v1, v2, Limits::Small);
 }
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.h	(revision 444)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.h	(revision 448)
@@ -20,5 +20,5 @@
 class VspKdTree;
 class AxisAlignedBox3;
-class VspBspLeaf;
+class BspLeaf;
 
 /**
@@ -57,7 +57,7 @@
 	*/
 	void  ComputeSampleContributions(const RayContainer &rays, 
-									 const bool castRays,
 									 int &sampleContributions,
-									 int &contributingSamples);
+									 int &contributingSamples,
+									 const bool castRays = true);
 
 
@@ -350,9 +350,9 @@
 	/** Merges view cells front and back leaf view cell.
 	*/
-	bool MergeVspBspLeafViewCells(VspBspLeaf *front, VspBspLeaf *back) const;
+	bool MergeVspBspLeafViewCells(BspLeaf *front, BspLeaf *back) const;
 
 	/** Returns true if front and back leaf should be merged.
 	*/
-	bool ShouldMerge(VspBspLeaf *front, VspBspLeaf *back) const;
+	bool ShouldMerge(BspLeaf *front, BspLeaf *back) const;
 
 	/// the view space partition BSP tree.
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp	(revision 445)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp	(revision 448)
@@ -13,4 +13,5 @@
 #include "Exporter.h"
 #include "Plane3.h"
+#include "ViewCellBsp.h"
 
 //-- static members
@@ -24,5 +25,4 @@
 const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
 
-int VspBspNode::sMailId = 1;
 
 int VspBspTree::sFrontId = 0; 
@@ -30,85 +30,1648 @@
 int VspBspTree::sFrontAndBackId = 0;
 
+
 /****************************************************************/
-/*                  class VspBspNode implementation                */
+/*                  class VspBspTree implementation             */
 /****************************************************************/
 
-VspBspNode::VspBspNode(): 
-mParent(NULL)
-{}
-
-VspBspNode::VspBspNode(VspBspInterior *parent): 
-mParent(parent)
-{}
-
-
-bool VspBspNode::IsRoot() const 
-{
-	return mParent == NULL;
-}
-
-VspBspInterior *VspBspNode::GetParent() 
-{ 
-	return mParent; 
-}
-
-void VspBspNode::SetParent(VspBspInterior *parent)
-{
-	mParent = parent;
-}
-
-/****************************************************************/
-/*              class VspBspInterior implementation                */
-/****************************************************************/
-
-
-VspBspInterior::VspBspInterior(const Plane3 &plane): 
-mPlane(plane), mFront(NULL), mBack(NULL)
-{}
-
-VspBspInterior::~VspBspInterior() 
-{ 
-	DEL_PTR(mFront); 
-	DEL_PTR(mBack);
-}
-
-bool VspBspInterior::IsLeaf() const
-{ 
-	return false; 
-}
-
-VspBspNode *VspBspInterior::GetBack() 
-{
-	return mBack;
-}
-
-VspBspNode *VspBspInterior::GetFront() 
-{
-	return mFront;
-}
-
-Plane3 *VspBspInterior::GetPlane()
-{
-	return &mPlane;
-}
-
-void VspBspInterior::ReplaceChildLink(VspBspNode *oldChild, VspBspNode *newChild) 
-{
-	if (mBack == oldChild)
-		mBack = newChild;
+VspBspTree::VspBspTree(): 
+mRoot(NULL),
+mPvsUseArea(true)
+{
+	mRootCell = new BspViewCell();
+
+	Randomize(); // initialise random generator for heuristics
+
+	//-- termination criteria for autopartition
+	environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
+	environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
+	environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
+	environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);	
+	environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
+	environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
+
+	//-- factors for bsp tree split plane heuristics
+	environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
+	environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
+	environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
+
+	//-- termination criteria for axis aligned split
+	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAaCtDivCi);
+	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
+	
+	environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays", 
+							 mTermMinRaysForAxisAligned);
+	
+	//-- partition criteria
+	environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
+	environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
+	environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
+	environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mSplitBorder);
+
+	environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
+	environment->GetIntValue("VspBspTree.maxTests", mMaxTests);
+
+    Debug << "BSP max depth: " << mTermMaxDepth << endl;
+	Debug << "BSP min PVS: " << mTermMinPvs << endl;
+	Debug << "BSP min area: " << mTermMinArea << endl;
+	Debug << "BSP min rays: " << mTermMinRays << endl;
+	Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
+	Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;
+
+	Debug << "Split plane strategy: ";
+	if (mSplitPlaneStrategy & RANDOM_POLYGON)
+		Debug << "random polygon ";
+	if (mSplitPlaneStrategy & AXIS_ALIGNED)
+		Debug << "axis aligned ";
+	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+		Debug << "least ray splits ";
+	if (mSplitPlaneStrategy & BALANCED_RAYS)
+		Debug << "balanced rays ";
+	if (mSplitPlaneStrategy & PVS)
+		Debug << "pvs";
+
+	Debug << endl;
+}
+
+
+const BspTreeStatistics &VspBspTree::GetStatistics() const 
+{
+	return mStat;
+}
+
+
+VspBspTree::~VspBspTree()
+{
+	DEL_PTR(mRoot);
+	DEL_PTR(mRootCell);
+}
+
+int VspBspTree::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(mEpsilon))
+		{
+			poly->mParent = parent; // set parent intersectable
+			polys.push_back(poly);
+		}
+		else
+			DEL_PTR(poly);
+	}
+	return (int)mesh->mFaces.size();
+}
+
+int VspBspTree::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 VspBspTree::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 VspBspTree::Construct(const VssRayContainer &sampleRays)
+{
+    mStat.nodes = 1;
+	mBox.Initialize(); 	// initialise BSP tree bounding box
+	
+	PolygonContainer polys;
+	RayInfoContainer *rays = new RayInfoContainer();
+
+	VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
+
+	long startTime = GetTime();
+
+	Debug << "**** Extracting polygons from rays ****\n";
+
+	Intersectable::NewMail();
+
+	//-- extract polygons intersected by the rays
+	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
+	{
+		VssRay *ray = *rit;
+	
+		if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
+		{
+			ray->mTerminationObject->Mail();
+			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
+			AddMeshToPolygons(obj->GetMesh(), polys, obj);
+		}
+
+		if (ray->mOriginObject && !ray->mOriginObject->Mailed())
+		{
+			ray->mOriginObject->Mail();
+			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
+			AddMeshToPolygons(obj->GetMesh(), polys, obj);
+		}
+	}
+
+	// compute bounding box
+	Polygon3::IncludeInBox(polys, mBox);
+
+	//-- store rays
+	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
+	{
+		VssRay *ray = *rit;
+		
+		float minT, maxT;
+
+		// TODO: not very efficient to implictly cast between rays types ...
+		if (mBox.GetRaySegment(*ray, minT, maxT))
+		{
+			float len = ray->Length();
+			
+			if (!len) 
+				len = Limits::Small;
+			
+			rays->push_back(RayInfo(ray, minT / len, maxT / len));
+		}
+	}
+
+	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);
+
+	// clean up polygons
+	CLEAR_CONTAINER(polys);
+}
+
+void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
+{
+	std::stack<VspBspTraversalData> tStack;
+
+	mRoot = new BspLeaf();
+
+	// constrruct root node geometry
+	BspNodeGeometry *geom = new BspNodeGeometry();
+	ConstructGeometry(mRoot, *geom);
+
+	VspBspTraversalData tData(mRoot, 
+							  new PolygonContainer(polys), 
+							  0,
+							  rays, 
+                              ComputePvsSize(*rays), 
+							  geom->GetArea(), 
+							  geom);
+
+	tStack.push(tData);
+
+	mStat.Start();
+	cout << "Contructing vsp bsp tree ... ";
+
+	long startTime = GetTime();
+	
+	while (!tStack.empty()) 
+	{
+		tData = tStack.top();
+
+	    tStack.pop();
+
+		// subdivide leaf node
+		BspNode *r = Subdivide(tStack, tData);
+
+		if (r == mRoot)
+			Debug << "BSP tree construction time spent at root: " 
+				  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+	}
+
+	cout << "finished\n";
+
+	mStat.Stop();
+}
+
+bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
+{
+	return 
+		(((int)data.mRays->size() <= mTermMinRays) ||
+		 (data.mPvs <= mTermMinPvs) ||
+		 (data.mArea <= mTermMinArea) ||
+		// (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
+		 (data.mDepth >= mTermMaxDepth));
+}
+
+BspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, 
+							   VspBspTraversalData &tData)
+{
+	//-- terminate traversal  
+	if (TerminationCriteriaMet(tData))		
+	{
+		BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
+	
+		BspViewCell *viewCell = new BspViewCell();
+		
+		leaf->SetViewCell(viewCell);
+		viewCell->mBspLeaves.push_back(leaf);
+
+		//-- add pvs
+		if (viewCell != mRootCell)
+		{
+			int conSamp = 0, sampCon = 0;
+			AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
+			
+			mStat.contributingSamples += conSamp;
+			mStat.sampleContributions += sampCon;
+		}
+
+		EvaluateLeafStats(tData);
+	
+		//-- clean up
+		
+		DEL_PTR(tData.mPolygons);
+		DEL_PTR(tData.mRays);
+		DEL_PTR(tData.mGeometry);
+
+		return leaf;
+	}
+
+	//-- continue subdivision
+	PolygonContainer coincident;
+	
+	VspBspTraversalData tFrontData(new PolygonContainer(), 
+								   tData.mDepth + 1,
+								   new RayInfoContainer(), 
+								   new BspNodeGeometry());
+
+	VspBspTraversalData tBackData(new PolygonContainer(), 
+						 		  tData.mDepth + 1,
+								  new RayInfoContainer(), 
+								  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
+
+	// 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;
+}
+
+BspInterior *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
+									   VspBspTraversalData &frontData,
+									   VspBspTraversalData &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->GetPlane(), 
+			  *tData.mRays, 
+			  *frontData.mRays, 
+			  *backData.mRays);
+	
+	// subdivide polygons with plane
+	mStat.splits += SplitPolygons(interior->GetPlane(),
+								  *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, 
+									   interior->GetPlane(),
+									   mBox,
+									   mEpsilon);
+	
+		
+		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->GetFront();
+	backData.mNode = interior->GetBack();
+	
+	//DEL_PTR(leaf);
+	return interior;
+}
+
+void VspBspTree::AddToPvs(BspLeaf *leaf,
+						  const RayInfoContainer &rays, 
+						  int &sampleContributions,
+						  int &contributingSamples)
+{
+	sampleContributions = 0;
+	contributingSamples = 0;
+
+    RayInfoContainer::const_iterator it, it_end = rays.end();
+
+	ViewCell *vc = leaf->GetViewCell();
+
+	// add contributions from samples to the PVS
+	for (it = rays.begin(); it != it_end; ++ it)
+	{
+		int contribution = 0;
+		VssRay *ray = (*it).mRay;
+			
+		if (ray->mTerminationObject)
+			contribution += vc->GetPvs().AddSample(ray->mTerminationObject);
+		
+		if (ray->mOriginObject)
+			contribution += vc->GetPvs().AddSample(ray->mOriginObject);
+
+		if (contribution)
+		{
+			sampleContributions += contribution;
+			++ contributingSamples;
+		}
+		//leaf->mVssRays.push_back(ray);
+	}
+}
+
+void VspBspTree::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 VspBspTree::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 = mAaCtDivCi + 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 VspBspTree::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 VspBspTree::SelectPlane(BspLeaf *leaf, VspBspTraversalData &data)
+{
+	if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
+		((int)data.mRays->size() > mTermMinRaysForAxisAligned))
+	{
+		Plane3 plane;
+		if (SelectAxisAlignedPlane(plane, *data.mPolygons)) // TODO: should be rays
+			return plane;
+	}
+
+	// simplest strategy: just take next polygon
+	if (mSplitPlaneStrategy & RANDOM_POLYGON)
+	{
+        if (!data.mPolygons->empty())
+		{
+			const int randIdx = Random((int)data.mPolygons->size());
+			Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
+			return nextPoly->GetSupportingPlane();
+		}
+		else
+		{
+			//-- choose plane on midpoint of a ray
+			const int candidateIdx = Random((int)data.mRays->size());
+									
+			const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
+			const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();
+
+			const Vector3 pt = (maxPt + minPt) * 0.5;
+
+			const Vector3 normal = (*data.mRays)[candidateIdx].mRay->GetDir();
+			
+			return Plane3(normal, pt);
+		}
+
+		return Plane3();
+	}
+
+	// use heuristics to find appropriate plane
+	return SelectPlaneHeuristics(leaf, data); 
+}
+
+Plane3 VspBspTree::SelectPlaneHeuristics(BspLeaf *leaf, VspBspTraversalData &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 different two methods chosen with
+	// equal probability
+	
+	// take 3 ray endpoints, where two are minimum and one a maximum
+	// point or the other way round
+	for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
+	{
+		candidateIdx = Random((int)data.mRays->size());
+	
+		RayInfo rayInf = (*data.mRays)[candidateIdx];
+
+		const Vector3 minPt = rayInf.ExtrapOrigin();
+		const Vector3 maxPt = rayInf.ExtrapTermination();
+
+		const Vector3 pt = (maxPt + minPt) * 0.5;
+		const Vector3 normal = Normalize(rayInf.mRay->GetDir());
+
+		plane = Plane3(normal, pt);
+
+		const float candidateCost = SplitPlaneCost(plane, data);
+
+		if (candidateCost < lowestCost)
+		{
+			bestPlane = plane;	
+			lowestCost = candidateCost;
+		}
+	}
+
+	// take plane normal as plane normal and the midpoint of the ray.
+	// PROBLEM: does not resemble any point where visibility is likely to change
+	//Debug << "lowest2: " << 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)data.mRays->size() * 2);
+				
+			if (idx[j] >= (int)data.mRays->size())
+			{
+				idx[j] -= (int)data.mRays->size();
+				
+				chooseMin = (cminT < 2);
+			}
+			else
+				chooseMin = (cmaxT < 2);
+
+			RayInfo rayInf = (*data.mRays)[idx[j]];
+			pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
+		}	
+			
+		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 VspBspTree::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());
+}
+
+
+inline void VspBspTree::GenerateUniqueIdsForPvs()
+{
+	Intersectable::NewMail(); sBackId = ViewCell::sMailId;
+	Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
+	Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
+}
+
+float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
+								 const VspBspTraversalData &data)
+{
+	float val = 0;
+
+	float sumBalancedRays = 0;
+	float sumRaySplits = 0;
+	
+	int frontPvs = 0;
+	int backPvs = 0;
+
+	// probability that view point lies in child
+	float pOverall = 0;
+	float pFront = 0;
+	float pBack = 0;
+
+	const bool pvsUseLen = false;
+
+	if (mSplitPlaneStrategy & PVS)
+	{
+		// create unique ids for pvs heuristics
+		GenerateUniqueIdsForPvs();
+	
+		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;
+		
+			data.mGeometry->SplitGeometry(frontCell, 
+										  backCell, 
+										  candidatePlane,
+										  mBox,
+										  mEpsilon);
+		
+			pFront = frontCell.GetArea();
+			pBack = backCell.GetArea();
+
+			pOverall = data.mArea;
+		}
+	}
+		
+	int limit;
+	bool useRand;
+
+	// choose test polyongs randomly if over threshold
+	if ((int)data.mRays->size() > mMaxTests)
+	{
+		useRand = true;
+		limit = mMaxTests;
+	}
 	else
-		mFront = newChild;
-}
-
-void VspBspInterior::SetupChildLinks(VspBspNode *b, VspBspNode *f) 
-{
-    mBack = b;
-    mFront = f;
-}
-
-int VspBspInterior::SplitPolygons(PolygonContainer &polys, 
-								  PolygonContainer &frontPolys, 
-								  PolygonContainer &backPolys, 
-								  PolygonContainer &coincident)
+	{
+		useRand = false;
+		limit = (int)data.mRays->size();
+	}
+
+	for (int i = 0; i < limit; ++ i)
+	{
+		const int testIdx = useRand ? Random(limit) : i;
+		RayInfo rayInf = (*data.mRays)[testIdx];
+
+		VssRay *ray = rayInf.mRay;
+		const int cf = rayInf.ComputeRayIntersection(candidatePlane, ray->mT);
+
+		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+		{
+			sumBalancedRays += cf;
+		}
+		
+		if (mSplitPlaneStrategy & BALANCED_RAYS)
+		{
+			if (cf == 0)
+				++ sumRaySplits;
+		}
+
+		if (mSplitPlaneStrategy & PVS)
+		{
+			// in case the ray intersects an object
+			// assure that we only count the object 
+			// once for the front and once for the back side of the plane
+			
+			// add the termination object
+			AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
+			
+			// add the source object
+			AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
+			
+			// use number or length of rays to approximate volume
+			if (!mPvsUseArea)
+			{
+				float len = 1;
+
+				if (pvsUseLen) // use length of rays
+					len = rayInf.SqrSegmentLength();
+			
+				pOverall += len;
+
+				if (cf == 1)
+					pFront += len;
+				if (cf == -1)
+					pBack += len;
+				if (cf == 0)
+				{
+					// use length of rays to approximate volume
+					if (pvsUseLen) 
+					{
+						float newLen = len * 
+							(rayInf.GetMaxT() - rayInf.mRay->mT) / 
+							(rayInf.GetMaxT() - rayInf.GetMinT());
+		
+						if (candidatePlane.Side(rayInf.ExtrapOrigin()) <= 0)
+						{
+							pBack += newLen;
+							pFront += len - newLen;
+						}
+						else
+						{
+							pFront += newLen;
+							pBack += len - newLen;
+						}
+					}
+					else
+					{
+						++ pFront;
+						++ pBack;
+					}
+				}
+			}
+		}
+	}
+
+	const float raysSize = (float)data.mRays->size() + Limits::Small;
+
+	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+		val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
+
+	if (mSplitPlaneStrategy & BALANCED_RAYS)
+		val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
+
+	const float denom = pOverall * (float)data.mPvs * 2.0f + Limits::Small;
+
+	if (mSplitPlaneStrategy & PVS)
+	{
+		val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
+
+		// give penalty to unbalanced split
+		if (1)
+		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 VspBspTree::AddObjToPvs(Intersectable *obj,
+							 const int cf,
+							 int &frontPvs,
+							 int &backPvs) const
+{
+	if (!obj)
+		return;
+	// TODO: does this really belong to no pvs?
+	//if (cf == Ray::COINCIDENT) return;
+
+	// object belongs to both PVS
+	if (cf >= 0)
+	{
+		if ((obj->mMailbox != sFrontId) && 
+			(obj->mMailbox != sFrontAndBackId))
+		{
+			++ frontPvs;
+
+			if (obj->mMailbox == sBackId)
+				obj->mMailbox = sFrontAndBackId;	
+			else
+				obj->mMailbox = sFrontId;								
+		}
+	}
+	
+	if (cf <= 0)
+	{
+		if ((obj->mMailbox != sBackId) &&
+			(obj->mMailbox != sFrontAndBackId))
+		{
+			++ backPvs;
+
+			if (obj->mMailbox == sFrontId)
+				obj->mMailbox = sFrontAndBackId; 
+			else
+				obj->mMailbox = sBackId;				
+		}
+	}
+}
+
+void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
+{
+	stack<BspNode *> nodeStack;
+	nodeStack.push(mRoot);
+ 
+	while (!nodeStack.empty()) 
+	{
+		BspNode *node = nodeStack.top();
+    
+		nodeStack.pop();
+    
+		if (node->IsLeaf()) 
+		{
+			BspLeaf *leaf = (BspLeaf *)node;		
+			leaves.push_back(leaf);
+		} 
+		else 
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+			nodeStack.push(interior->GetBack());
+			nodeStack.push(interior->GetFront());
+		}
+	}
+}
+
+AxisAlignedBox3 VspBspTree::GetBoundingBox() const
+{
+	return mBox;
+}
+
+BspNode *VspBspTree::GetRoot() const
+{
+	return mRoot;
+}
+
+void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
+{
+	// the node became a leaf -> evaluate stats for leafs
+	BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
+
+	// store maximal and minimal depth
+	if (data.mDepth > mStat.maxDepth)
+		mStat.maxDepth = data.mDepth;
+
+	if (data.mDepth < mStat.minDepth)
+		mStat.minDepth = data.mDepth;
+
+	// accumulate depth to compute average depth
+	mStat.accumDepth += data.mDepth;
+	
+
+	if (data.mDepth >= mTermMaxDepth)
+		++ mStat.maxDepthNodes;
+
+	if (data.mPvs < mTermMinPvs)
+		++ mStat.minPvsNodes;
+
+	if ((int)data.mRays->size() < mTermMinRays)
+		++ mStat.minRaysNodes;
+
+	if (data.GetAvgRayContribution() > mTermMaxRayContribution)
+		++ mStat.maxRayContribNodes;
+	
+	if (data.mGeometry->GetArea() <= mTermMinArea) 
+		++ mStat.minAreaNodes;
+
+#ifdef _DEBUG
+	Debug << "BSP stats: "
+		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
+		  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
+		  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
+		  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), " 
+		  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
+		  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
+#endif
+}
+
+int VspBspTree::CastRay(Ray &ray)
+{
+	int hits = 0;
+  
+	stack<BspRayTraversalData> tStack;
+  
+	float maxt, mint;
+
+	if (!mBox.GetRaySegment(ray, mint, maxt))
+		return 0;
+
+	Intersectable::NewMail();
+
+	Vector3 entp = ray.Extrap(mint);
+	Vector3 extp = ray.Extrap(maxt);
+  
+	BspNode *node = mRoot;
+	BspNode *farChild = NULL;
+	
+	while (1)
+	{
+		if (!node->IsLeaf()) 
+		{
+			BspInterior *in = dynamic_cast<BspInterior *>(node);
+
+			Plane3 splitPlane = in->GetPlane();
+			const int entSide = splitPlane.Side(entp);
+			const int extSide = splitPlane.Side(extp);
+
+			if (entSide < 0)
+			{
+				node = in->GetBack();
+
+				if(extSide <= 0) // plane does not split ray => no far child
+					continue;
+					
+				farChild = in->GetFront(); // plane splits ray
+
+			} else if (entSide > 0)
+			{
+				node = in->GetFront();
+
+				if (extSide >= 0) // plane does not split ray => no far child
+					continue;
+
+				farChild = in->GetBack(); // plane splits ray			
+			}
+			else // ray and plane are coincident
+			{
+				// WHAT TO DO IN THIS CASE ?
+				//break;
+				node = in->GetFront();
+				continue;
+			}
+
+			// push data for far child
+			tStack.push(BspRayTraversalData(farChild, extp, maxt));
+
+			// find intersection of ray segment with plane
+			float t;
+			extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
+			maxt *= t;
+			
+		} else // reached leaf => intersection with view cell
+		{
+			BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
+      
+			if (!leaf->GetViewCell()->Mailed())
+			{
+				//ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
+				leaf->GetViewCell()->Mail();
+				++ hits;
+			}
+			
+   			//-- fetch the next far child from the stack
+			if (tStack.empty())
+				break;
+      
+			entp = extp;
+			mint = maxt; // NOTE: need this?
+
+			if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
+				break;
+
+			BspRayTraversalData &s = tStack.top();
+
+			node = s.mNode;
+			extp = s.mExitPoint;
+			maxt = s.mMaxT;
+
+			tStack.pop();
+		}
+	}
+
+	return hits;
+}
+
+bool VspBspTree::Export(const string filename)
+{
+	Exporter *exporter = Exporter::GetExporter(filename);
+
+	if (exporter) 
+	{
+		//exporter->ExportVspBspTree(*this);
+		return true;
+	}	
+
+	return false;
+}
+
+void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
+{
+	stack<BspNode *> nodeStack;
+	nodeStack.push(mRoot);
+
+	ViewCell::NewMail();
+
+	while (!nodeStack.empty()) 
+	{
+		BspNode *node = nodeStack.top();
+		nodeStack.pop();
+
+		if (node->IsLeaf()) 
+		{
+			ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
+
+			if (!viewCell->Mailed()) 
+			{
+				viewCell->Mail();
+				viewCells.push_back(viewCell);
+			}
+		}
+		else 
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+			nodeStack.push(interior->GetFront());
+			nodeStack.push(interior->GetBack());
+		}
+	}
+}
+
+void VspBspTree::EvaluateViewCellsStats(BspViewCellsStatistics &stat) const
+{
+	stat.Reset();
+
+	stack<BspNode *> nodeStack;
+	nodeStack.push(mRoot);
+
+	ViewCell::NewMail();
+
+	// exclude root cell
+	mRootCell->Mail();
+
+	while (!nodeStack.empty()) 
+	{
+		BspNode *node = nodeStack.top();
+		nodeStack.pop();
+
+		if (node->IsLeaf()) 
+		{
+			++ stat.bspLeaves;
+
+			BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
+
+			if (!viewCell->Mailed()) 
+			{
+				viewCell->Mail();
+				
+				++ stat.viewCells;
+				const int pvsSize = viewCell->GetPvs().GetSize();
+
+                stat.pvs += pvsSize;
+
+				if (pvsSize < 1)
+					++ stat.emptyPvs;
+
+				if (pvsSize > stat.maxPvs)
+					stat.maxPvs = pvsSize;
+
+				if (pvsSize < stat.minPvs)
+					stat.minPvs = pvsSize;
+
+				if ((int)viewCell->mBspLeaves.size() > stat.maxBspLeaves)
+					stat.maxBspLeaves = (int)viewCell->mBspLeaves.size();		
+			}
+		}
+		else 
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+			nodeStack.push(interior->GetFront());
+			nodeStack.push(interior->GetBack());
+		}
+	}
+}
+
+BspTreeStatistics &VspBspTree::GetStat()
+{
+	return mStat;
+}
+
+float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
+{
+	float len = 0;
+
+	RayInfoContainer::const_iterator it, it_end = rays.end();
+
+	for (it = rays.begin(); it != it_end; ++ it)
+		len += (*it).SegmentLength();
+
+	return len;
+}
+
+int VspBspTree::SplitRays(const Plane3 &plane,
+						  RayInfoContainer &rays, 
+						  RayInfoContainer &frontRays, 
+						  RayInfoContainer &backRays)
+{
+	int splits = 0;
+
+	while (!rays.empty())
+	{
+		RayInfo bRay = rays.back();
+		VssRay *ray = bRay.mRay;
+
+		rays.pop_back();
+	
+		// get classification and receive new t
+		const int cf = bRay.ComputeRayIntersection(plane, ray->mT);
+
+		switch (cf)
+		{
+		case -1:
+			backRays.push_back(bRay);
+			break;
+		case 1:
+			frontRays.push_back(bRay);
+			break;
+		case 0: 
+			//-- split ray
+
+			// if start point behind plane
+			if (plane.Side(bRay.ExtrapOrigin()) <= 0)
+			{	
+				backRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
+				frontRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
+			}
+			else
+			{
+				frontRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
+				backRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
+			}
+			break;
+		default:
+			Debug << "Should not come here 4" << endl;
+			break;
+		}
+	}
+
+	return splits;
+}
+
+void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
+{
+	BspNode *lastNode;
+
+	do
+	{
+		lastNode = n;
+
+		// want to get planes defining geometry of this node => don't take
+		// split plane of node itself
+		n = n->GetParent();
+		
+		if (n)
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(n);
+			Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
+
+            if (interior->GetFront() != lastNode)
+				halfSpace.ReverseOrientation();
+
+			halfSpaces.push_back(halfSpace);
+		}
+	}
+	while (n);
+}
+
+void VspBspTree::ConstructGeometry(BspNode *n, 
+								   BspNodeGeometry &cell) const
+{
+	PolygonContainer polys;
+	ConstructGeometry(n, polys);
+	cell.mPolys = polys;
+}
+
+void VspBspTree::ConstructGeometry(BspNode *n, 
+								   PolygonContainer &cell) const
+{
+	vector<Plane3> halfSpaces;
+	ExtractHalfSpaces(n, halfSpaces);
+
+	PolygonContainer candidatePolys;
+
+	// bounded planes are added to the polygons (reverse polygons
+	// as they have to be outfacing
+	for (int i = 0; i < (int)halfSpaces.size(); ++ i)
+	{
+		Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
+		
+		if (p->Valid(mEpsilon))
+		{
+			candidatePolys.push_back(p->CreateReversePolygon());
+			DEL_PTR(p);
+		}
+	}
+
+	// add faces of bounding box (also could be faces of the cell)
+	for (int i = 0; i < 6; ++ i)
+	{
+		VertexContainer vertices;
+	
+		for (int j = 0; j < 4; ++ j)
+			vertices.push_back(mBox.GetFace(i).mVertices[j]);
+
+		candidatePolys.push_back(new Polygon3(vertices));
+	}
+
+	for (int i = 0; i < (int)candidatePolys.size(); ++ i)
+	{
+		// polygon is split by all other planes
+		for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
+		{
+			if (i == j) // polygon and plane are coincident
+				continue;
+
+			VertexContainer splitPts;
+			Polygon3 *frontPoly, *backPoly;
+
+			const int cf = 
+				candidatePolys[i]->ClassifyPlane(halfSpaces[j],
+												 mEpsilon);
+			
+			switch (cf)
+			{
+				case Polygon3::SPLIT:
+					frontPoly = new Polygon3();
+					backPoly = new Polygon3();
+
+					candidatePolys[i]->Split(halfSpaces[j], 
+											 *frontPoly, 
+											 *backPoly,
+											 mEpsilon);
+
+					DEL_PTR(candidatePolys[i]);
+
+					if (frontPoly->Valid(mEpsilon))
+						candidatePolys[i] = frontPoly;
+					else
+						DEL_PTR(frontPoly);
+
+					DEL_PTR(backPoly);
+					break;
+				case Polygon3::BACK_SIDE:
+					DEL_PTR(candidatePolys[i]);
+					break;
+				// just take polygon as it is
+				case Polygon3::FRONT_SIDE:
+				case Polygon3::COINCIDENT:
+				default:
+					break;
+			}
+		}
+		
+		if (candidatePolys[i])
+			cell.push_back(candidatePolys[i]);
+	}
+}
+
+void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &vcGeom) const
+{
+	vector<BspLeaf *> leaves = vc->mBspLeaves;
+
+	vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
+
+	for (it = leaves.begin(); it != it_end; ++ it)
+		ConstructGeometry(*it, vcGeom);
+}
+
+int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors, 
+						   const bool onlyUnmailed) const
+{
+	PolygonContainer cell;
+
+	ConstructGeometry(n, cell);
+
+	stack<BspNode *> nodeStack;
+	nodeStack.push(mRoot);
+		
+	// planes needed to verify that we found neighbor leaf.
+	vector<Plane3> halfSpaces;
+	ExtractHalfSpaces(n, halfSpaces);
+
+	while (!nodeStack.empty()) 
+	{
+		BspNode *node = nodeStack.top();
+		nodeStack.pop();
+
+		if (node->IsLeaf())
+		{
+            if (node != n && (!onlyUnmailed || !node->Mailed()))
+			{
+				// test all planes of current node if candidate really
+				// is neighbour
+				PolygonContainer neighborCandidate;
+				ConstructGeometry(node, neighborCandidate);
+				
+				bool isAdjacent = true;
+				for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
+				{
+					const int cf = 
+						Polygon3::ClassifyPlane(neighborCandidate, 
+												halfSpaces[i],
+												mEpsilon);
+
+					if (cf == Polygon3::BACK_SIDE)
+						isAdjacent = false;
+				}
+
+				if (isAdjacent)
+					neighbors.push_back(dynamic_cast<BspLeaf *>(node));
+
+				CLEAR_CONTAINER(neighborCandidate);
+			}
+		}
+		else 
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+	
+			const int cf = Polygon3::ClassifyPlane(cell, 
+												   interior->GetPlane(), 
+												   mEpsilon);
+
+			if (cf == Polygon3::FRONT_SIDE)
+				nodeStack.push(interior->GetFront());
+			else
+				if (cf == Polygon3::BACK_SIDE)
+					nodeStack.push(interior->GetBack());
+				else 
+				{
+					// random decision
+					nodeStack.push(interior->GetBack());
+					nodeStack.push(interior->GetFront());
+				}
+		}
+	}
+	
+	CLEAR_CONTAINER(cell);
+	return (int)neighbors.size();
+}
+
+BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
+{
+    stack<BspNode *> nodeStack;
+  	nodeStack.push(mRoot);
+	
+	int mask = rand();
+  
+	while (!nodeStack.empty()) 
+	{
+		BspNode *node = nodeStack.top();
+		nodeStack.pop();
+	  
+		if (node->IsLeaf()) 
+		{
+			return dynamic_cast<BspLeaf *>(node);
+		} 
+		else 
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+			
+			BspNode *next;
+	
+			PolygonContainer cell;
+
+			// todo: not very efficient: constructs full cell everytime
+			ConstructGeometry(interior, cell);
+
+			const int cf = Polygon3::ClassifyPlane(cell, halfspace, mEpsilon);
+
+			if (cf == Polygon3::BACK_SIDE)
+				next = interior->GetFront();
+			else
+				if (cf == Polygon3::FRONT_SIDE)
+					next = interior->GetFront();
+			else 
+			{
+				// random decision
+				if (mask & 1)
+					next = interior->GetBack();
+				else
+					next = interior->GetFront();
+				mask = mask >> 1;
+			}
+
+			nodeStack.push(next);
+		}
+	}
+	
+	return NULL;
+}
+
+BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
+{
+	stack<BspNode *> nodeStack;
+	
+	nodeStack.push(mRoot);
+
+	int mask = rand();
+	
+	while (!nodeStack.empty()) 
+	{
+		BspNode *node = nodeStack.top();
+		nodeStack.pop();
+		
+		if (node->IsLeaf()) 
+		{
+			if ( (!onlyUnmailed || !node->Mailed()) )
+				return dynamic_cast<BspLeaf *>(node);
+		}
+		else 
+		{
+			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+			// random decision
+			if (mask & 1)
+				nodeStack.push(interior->GetBack());
+			else
+				nodeStack.push(interior->GetFront());
+
+			mask = mask >> 1;
+		}
+	}
+	
+	return NULL;
+}
+
+int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
+{
+	int pvsSize = 0;
+
+	RayInfoContainer::const_iterator rit, rit_end = rays.end();
+
+	Intersectable::NewMail();
+
+	for (rit = rays.begin(); rit != rays.end(); ++ rit)
+	{
+		VssRay *ray = (*rit).mRay;
+		
+		if (ray->mOriginObject)
+		{
+			if (!ray->mOriginObject->Mailed())
+			{
+				ray->mOriginObject->Mail();
+				++ pvsSize;
+			}
+		}
+		if (ray->mTerminationObject)
+		{
+			if (!ray->mTerminationObject->Mailed())
+			{
+				ray->mTerminationObject->Mail();
+				++ pvsSize;
+			}
+		}
+	}
+
+	return pvsSize;
+}
+
+float VspBspTree::GetEpsilon() const
+{
+	return mEpsilon;
+}
+
+BspViewCell *VspBspTree::GetRootCell() const
+{
+	return mRootCell;
+}
+
+int VspBspTree::SplitPolygons(const Plane3 &plane,
+							  PolygonContainer &polys, 
+							  PolygonContainer &frontPolys, 
+							  PolygonContainer &backPolys, 
+							  PolygonContainer &coincident) const
 {
 	int splits = 0;
@@ -120,5 +1683,5 @@
 
 		// classify polygon
-		const int cf = poly->ClassifyPlane(mPlane);
+		const int cf = poly->ClassifyPlane(plane, mEpsilon);
 
 		switch (cf)
@@ -146,1822 +1709,2 @@
 	return splits;
 }
-
-/****************************************************************/
-/*                  class VspBspLeaf implementation                */
-/****************************************************************/
-VspBspLeaf::VspBspLeaf(): mViewCell(NULL)
-{
-}
-
-VspBspLeaf::VspBspLeaf(BspViewCell *viewCell): 
-mViewCell(viewCell)
-{
-}
-
-VspBspLeaf::VspBspLeaf(VspBspInterior *parent): 
-VspBspNode(parent), mViewCell(NULL)
-{}
-
-VspBspLeaf::VspBspLeaf(VspBspInterior *parent, BspViewCell *viewCell): 
-VspBspNode(parent), mViewCell(viewCell)
-{
-}
-
-BspViewCell *VspBspLeaf::GetViewCell() const
-{
-	return mViewCell;
-}
-
-void VspBspLeaf::SetViewCell(BspViewCell *viewCell)
-{
-	mViewCell = viewCell;
-}
-
-bool VspBspLeaf::IsLeaf() const 
-{ 
-	return true; 
-}
-
-void VspBspLeaf::AddToPvs(const RayInfoContainer &rays, 
-						  int &sampleContributions,
-						  int &contributingSamples,
-						  bool storeLeavesWithRays)
-{
-	sampleContributions = 0;
-	contributingSamples = 0;
-
-    RayInfoContainer::const_iterator it, it_end = rays.end();
-
-	// add contributions from samples to the PVS
-	for (it = rays.begin(); it != it_end; ++ it)
-	{
-		int contribution = 0;
-		VssRay *ray = (*it).mRay;
-			
-		if (ray->mTerminationObject)
-			contribution += mViewCell->GetPvs().AddSample(ray->mTerminationObject);
-		
-		if (ray->mOriginObject)
-			contribution += mViewCell->GetPvs().AddSample(ray->mOriginObject);
-
-		if (contribution > 0)
-		{
-			sampleContributions += contribution;
-			++ contributingSamples;
-		}
-	}
-}
-
-
-/****************************************************************/
-/*                  class VspBspTree implementation             */
-/****************************************************************/
-
-VspBspTree::VspBspTree(): 
-mRoot(NULL),
-mPvsUseArea(true)
-{
-	mRootCell = new BspViewCell();
-
-	Randomize(); // initialise random generator for heuristics
-
-	//-- termination criteria for autopartition
-	environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
-	environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
-	environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
-	environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);	
-	environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
-	environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
-
-	//-- factors for bsp tree split plane heuristics
-	environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
-	environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
-	environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
-
-	//-- termination criteria for axis aligned split
-	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAaCtDivCi);
-	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
-	
-	environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays", 
-							 mTermMinRaysForAxisAligned);
-	
-	//-- partition criteria
-	environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
-	environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
-	environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
-	environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mSplitBorder);
-
-	environment->GetFloatValue("VspBspTree.Construction.sideTolerance", Vector3::sDistTolerance);
-	Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
-
-    Debug << "BSP max depth: " << mTermMaxDepth << endl;
-	Debug << "BSP min PVS: " << mTermMinPvs << endl;
-	Debug << "BSP min area: " << mTermMinArea << endl;
-	Debug << "BSP min rays: " << mTermMinRays << endl;
-	Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
-	Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;
-
-	Debug << "Split plane strategy: ";
-	if (mSplitPlaneStrategy & RANDOM_POLYGON)
-		Debug << "random polygon ";
-	if (mSplitPlaneStrategy & AXIS_ALIGNED)
-		Debug << "axis aligned ";
-	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		Debug << "least ray splits ";
-	if (mSplitPlaneStrategy & BALANCED_RAYS)
-		Debug << "balanced rays ";
-	if (mSplitPlaneStrategy & PVS)
-		Debug << "pvs";
-
-	Debug << endl;
-}
-
-
-const VspBspTreeStatistics &VspBspTree::GetStatistics() const 
-{
-	return mStat;
-}
-
-void VspBspTreeStatistics::Print(ostream &app) const
-{
-	app << "===== VspBspTree 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 maximum depth )\n" 
-		<<	maxDepthNodes * 100 / (double)Leaves() << endl;
-
-	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
-		<<	maxDepthNodes * 100 / (double)Leaves() << endl;
-
-	app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n" 
-		<< minPvsNodes * 100 / (double)Leaves() << endl;
-
-	app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n" 
-		<<	minRaysNodes * 100 / (double)Leaves() << endl;
-
-	app << "#N_PMINAREALEAVES  ( Percentage of leaves with mininum area )\n"
-		<< minAreaNodes * 100 / (double)Leaves() << endl;
-
-	app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n" 
-		<<	maxRayContribNodes * 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 VspBspTree statistics ==========\n";
-}
-
-
-VspBspTree::~VspBspTree()
-{
-	DEL_PTR(mRoot);
-	DEL_PTR(mRootCell);
-}
-
-int VspBspTree::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 VspBspTree::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 VspBspTree::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 VspBspTree::Construct(const VssRayContainer &sampleRays)
-{
-    mStat.nodes = 1;
-	mBox.Initialize(); 	// initialise BSP tree bounding box
-	
-	PolygonContainer polys;
-	RayInfoContainer *rays = new RayInfoContainer();
-
-	VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
-
-	long startTime = GetTime();
-
-	Debug << "**** Extracting polygons from rays ****\n";
-
-	//-- extract polygons intersected by the rays
-	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
-	{
-		VssRay *ray = *rit;
-	
-		if (ray->mTerminationObject)
-		{
-			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
-			AddMeshToPolygons(obj->GetMesh(), polys, obj);
-		}
-
-		if (ray->mOriginObject)
-		{
-			MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
-			AddMeshToPolygons(obj->GetMesh(), polys, obj);
-		}
-	}
-
-	//-- store rays
-	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
-	{
-		rays->push_back(RayInfo(*rit));
-	}
-
-	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(new PolygonContainer(polys), rays);
-
-	// clean up polygons because they are not deleted in the tree
-	CLEAR_CONTAINER(polys);
-}
-
-void VspBspTree::Construct(PolygonContainer *polys, RayInfoContainer *rays)
-{
-	std::stack<VspBspTraversalData> tStack;
-
-	mRoot = new VspBspLeaf();
-
-	VspBspNodeGeometry *cell = new VspBspNodeGeometry();
-	ConstructGeometry(mRoot, *cell);
-
-	VspBspTraversalData tData(mRoot, polys, 0, mRootCell, rays, 
-						   ComputePvsSize(*rays), cell->GetArea(), cell);
-
-	tStack.push(tData);
-
-	mStat.Start();
-	cout << "Contructing bsp tree ... ";
-	long startTime = GetTime();
-	while (!tStack.empty()) 
-	{
-		tData = tStack.top();
-
-	    tStack.pop();
-
-		// subdivide leaf node
-		VspBspNode *r = Subdivide(tStack, tData);
-
-		if (r == mRoot)
-			Debug << "BSP tree construction time spent at root: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
-	}
-
-	cout << "finished\n";
-
-	mStat.Stop();
-}
-
-bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
-{
-	return 
-		(((int)data.mRays->size() <= mTermMinRays) ||
-		 (data.mPvs <= mTermMinPvs) ||
-		 (data.mArea <= mTermMinArea) ||
-		 (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
-		 (data.mDepth >= mTermMaxDepth));
-}
-
-VspBspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, VspBspTraversalData &tData)
-{
-	//-- terminate traversal  
-	if (TerminationCriteriaMet(tData))		
-	{
-		VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(tData.mNode);
-	
-		BspViewCell *viewCell = new BspViewCell();
-		
-		leaf->SetViewCell(viewCell);
-		//TODO viewCell->mVspBspLeaves.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
-
-		DEL_PTR(tData.mPolygons);
-		DEL_PTR(tData.mRays);
-		DEL_PTR(tData.mGeometry);
-		
-		return leaf;
-	}
-
-	//-- continue subdivision
-	PolygonContainer coincident;
-	
-	VspBspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
-								new RayInfoContainer(), 0, 0, new VspBspNodeGeometry());
-	VspBspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
-							   new RayInfoContainer(), 0, 0, new VspBspNodeGeometry());
-
-	// create new interior node and two leaf nodes
-	VspBspInterior *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
-
-	// 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;
-}
-
-VspBspInterior *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
-									VspBspTraversalData &frontData,
-									VspBspTraversalData &backData,
-									PolygonContainer &coincident)
-{
-	mStat.nodes += 2;
-	
-	VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(tData.mNode);
-	// select subdivision plane
-	VspBspInterior *interior = 
-		new VspBspInterior(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
-
-	VspBspInterior *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 VspBspLeaf(interior), new VspBspLeaf(interior));
-	
-	frontData.mNode = interior->mFront;
-	backData.mNode = interior->mBack;
-	
-	//DEL_PTR(leaf);
-	return interior;
-}
-
-void VspBspTree::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 VspBspTree::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 = mAaCtDivCi + 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 VspBspTree::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 VspBspTree::SelectPlane(VspBspLeaf *leaf, VspBspTraversalData &data)
-{
-	if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
-		((int)data.mRays->size() > mTermMinRaysForAxisAligned))
-	{
-		Plane3 plane;
-		if (SelectAxisAlignedPlane(plane, *data.mPolygons)) // TODO: should be rays
-			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
-		{
-			//-- choose plane on midpoint of a ray
-			const int candidateIdx = Random((int)data.mRays->size());
-								
-			const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
-			const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();
-
-			const Vector3 pt = (maxPt + minPt) * 0.5;
-
-			const Vector3 normal =  (*data.mRays)[candidateIdx].mRay->GetDir();
-			
-			return Plane3(normal, pt);
-		}
-
-		return Plane3();
-	}
-
-	// use heuristics to find appropriate plane
-	return SelectPlaneHeuristics(leaf, data); 
-}
-
-Plane3 VspBspTree::SelectPlaneHeuristics(VspBspLeaf *leaf, VspBspTraversalData &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 different two methods chosen with
-	// equal probability
-	
-	RayInfoContainer *rays = data.mRays;
-	// take 3 ray endpoints, where two are minimum and one a maximum
-	// point or the other way round
-	for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
-	{
-		candidateIdx = Random((int)rays->size());
-	
-		RayInfo rayInf = (*rays)[candidateIdx];
-		const Vector3 minPt = rayInf.ExtrapOrigin();
-		const Vector3 maxPt = rayInf.ExtrapTermination();
-
-		const Vector3 pt = (maxPt + minPt) * 0.5;
-		const Vector3 normal = Normalize(rayInf.mRay->GetDir());
-
-		plane = Plane3(normal, pt);
-	
-		const float candidateCost = SplitPlaneCost(plane, data);
-
-		if (candidateCost < lowestCost)
-		{
-			bestPlane = plane;	
-			lowestCost = candidateCost;
-		}
-	}
-
-	// take plane normal as plane normal and the midpoint of the ray.
-	// PROBLEM: does not resemble any point where visibility is likely to change
-	//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);
-
-			RayInfo rayInf = (*rays)[idx[j]];
-			pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
-		}	
-			
-		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 VspBspTree::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());
-}
-
-
-bool VspBspTree::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;
-}
-
-inline void VspBspTree::GenerateUniqueIdsForPvs()
-{
-	Intersectable::NewMail(); sBackId = ViewCell::sMailId;
-	Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
-	Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
-}
-
-float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
-								 const VspBspTraversalData &data)
-{
-	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;
-
-	const bool pvsUseLen = false;
-
-	if (mSplitPlaneStrategy & PVS)
-	{
-		// create unique ids for pvs heuristics
-		GenerateUniqueIdsForPvs();
-	
-		if (mPvsUseArea) // use front and back cell areas to approximate volume
-		{	
-			// construct child geometry with regard to the candidate split plane
-			VspBspNodeGeometry frontCell;
-			VspBspNodeGeometry backCell;
-
-			data.mGeometry->SplitGeometry(frontCell, backCell, *this, candidatePlane);
-		
-			pFront = frontCell.GetArea();
-			pBack = backCell.GetArea();
-
-			pOverall = data.mArea;
-		}
-	}
-			
-	RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
-
-	for (rit = data.mRays->begin(); rit != data.mRays->end(); ++ rit)
-	{
-		VssRay *ray = (*rit).mRay;
-		const int cf = (*rit).ComputeRayIntersection(candidatePlane, (*rit).mRay->mT);
-
-		if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		{
-			sumBalancedRays += cf;
-		}
-		
-		if (mSplitPlaneStrategy & BALANCED_RAYS)
-		{
-			if (cf == 0)
-				++ sumRaySplits;
-		}
-
-		if (mSplitPlaneStrategy & PVS)
-		{
-			// in case the ray intersects an object
-			// assure that we only count the object 
-			// once for the front and once for the back side of the plane
-			
-			// add the termination object
-			AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
-			
-			// add the source object
-			AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
-			
-			// use number or length of rays to approximate volume
-			if (!mPvsUseArea)
-			{
-				float len = 1;
-
-				if (pvsUseLen) // use length of rays
-					len = (*rit).SqrSegmentLength();
-			
-				pOverall += len;
-
-				if (cf == 1)
-					pFront += len;
-				if (cf == -1)
-					pBack += len;
-				if (cf == 0)
-				{
-					// use length of rays to approximate volume
-					if (pvsUseLen) 
-					{
-						float newLen = len * 
-							((*rit).GetMaxT() - (*rit).mRay->mT) / 
-							((*rit).GetMaxT() - (*rit).GetMinT());
-		
-						if (candidatePlane.Side((*rit).ExtrapOrigin()) <= 0)
-						{
-							pBack += newLen;
-							pFront += len - newLen;
-						}
-						else
-						{
-							pFront += newLen;
-							pBack += len - newLen;
-						}
-					}
-					else
-					{
-						++ pFront;
-						++ pBack;
-					}
-				}
-			}
-		}
-	}
-
-	const float raysSize = (float)data.mRays->size() + Limits::Small;
-
-	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-		val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
-
-	if (mSplitPlaneStrategy & BALANCED_RAYS)
-		val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
-
-	const float denom = pOverall * (float)data.mPvs * 2.0f + Limits::Small;
-
-	if (mSplitPlaneStrategy & 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 VspBspTree::AddObjToPvs(Intersectable *obj,
-						  const int cf,
-						  int &frontPvs,
-						  int &backPvs) const
-{
-	if (!obj)
-		return;
-	// TODO: does this really belong to no pvs?
-	//if (cf == Ray::COINCIDENT) return;
-
-	// object belongs to both PVS
-	if (cf >= 0)
-	{
-		if ((obj->mMailbox != sFrontId) && 
-			(obj->mMailbox != sFrontAndBackId))
-		{
-			++ frontPvs;
-
-			if (obj->mMailbox == sBackId)
-				obj->mMailbox = sFrontAndBackId;	
-			else
-				obj->mMailbox = sFrontId;								
-		}
-	}
-	
-	if (cf <= 0)
-	{
-		if ((obj->mMailbox != sBackId) &&
-			(obj->mMailbox != sFrontAndBackId))
-		{
-			++ backPvs;
-
-			if (obj->mMailbox == sFrontId)
-				obj->mMailbox = sFrontAndBackId; 
-			else
-				obj->mMailbox = sBackId;				
-		}
-	}
-}
-
-void VspBspTree::CollectLeaves(vector<VspBspLeaf *> &leaves) const
-{
-	stack<VspBspNode *> nodeStack;
-	nodeStack.push(mRoot);
- 
-	while (!nodeStack.empty()) 
-	{
-		VspBspNode *node = nodeStack.top();
-    
-		nodeStack.pop();
-    
-		if (node->IsLeaf()) 
-		{
-			VspBspLeaf *leaf = (VspBspLeaf *)node;		
-			leaves.push_back(leaf);
-		} 
-		else 
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-			nodeStack.push(interior->GetBack());
-			nodeStack.push(interior->GetFront());
-		}
-	}
-}
-
-AxisAlignedBox3 VspBspTree::GetBoundingBox() const
-{
-	return mBox;
-}
-
-VspBspNode *VspBspTree::GetRoot() const
-{
-	return mRoot;
-}
-
-void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
-{
-	// the node became a leaf -> evaluate stats for leafs
-	VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(data.mNode);
-
-	// store maximal and minimal depth
-	if (data.mDepth > mStat.maxDepth)
-		mStat.maxDepth = data.mDepth;
-
-	if (data.mDepth < mStat.minDepth)
-		mStat.minDepth = data.mDepth;
-
-	// accumulate depth to compute average depth
-	mStat.accumDepth += data.mDepth;
-	
-
-	if (data.mDepth >= mTermMaxDepth)
-		++ mStat.maxDepthNodes;
-
-	if (data.mPvs < mTermMinPvs)
-		++ mStat.minPvsNodes;
-
-	if ((int)data.mRays->size() < mTermMinRays)
-		++ mStat.minRaysNodes;
-
-	if (data.GetAvgRayContribution() > mTermMaxRayContribution)
-		++ mStat.maxRayContribNodes;
-	
-	if (data.mGeometry->GetArea() <= mTermMinArea) 
-		++ mStat.minAreaNodes;
-
-#ifdef _DEBUG
-	Debug << "BSP stats: "
-		  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
-		  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
-		  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
-		  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << mTermMinPolys << "), "
-		  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), " 
-		  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
-		  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
-#endif
-}
-
-int VspBspTree::CastRay(Ray &ray)
-{
-	int hits = 0;
-  
-	stack<VspBspRayTraversalData> tStack;
-  
-	float maxt, mint;
-
-	if (!BoundRay(ray, mint, maxt))
-		return 0;
-
-	Intersectable::NewMail();
-
-	Vector3 entp = ray.Extrap(mint);
-	Vector3 extp = ray.Extrap(maxt);
-  
-	VspBspNode *node = mRoot;
-	VspBspNode *farChild = NULL;
-	
-	while (1)
-	{
-		if (!node->IsLeaf()) 
-		{
-			VspBspInterior *in = (VspBspInterior *) node;
-			
-			Plane3 *splitPlane = in->GetPlane();
-
-			int entSide = splitPlane->Side(entp);
-			int extSide = splitPlane->Side(extp);
-
-			Vector3 intersection;
-
-			if (entSide < 0)
-			{
-				node = in->GetBack();
-
-				if(extSide <= 0) // plane does not split ray => no far child
-					continue;
-					
-				farChild = in->GetFront(); // plane splits ray
-
-			} else if (entSide > 0)
-			{
-				node = in->GetFront();
-
-				if (extSide >= 0) // plane does not split ray => no far child
-					continue;
-
-				farChild = in->GetBack(); // plane splits ray			
-			}
-			else // ray and plane are coincident
-			{
-				// WHAT TO DO IN THIS CASE ?
-				//break;
-				node = in->GetFront();
-				continue;
-			}
-
-			// push data for far child
-			tStack.push(VspBspRayTraversalData(farChild, extp, maxt));
-
-			// find intersection of ray segment with plane
-			float t;
-			extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &t);
-			maxt *= t;
-			
-		} else // reached leaf => intersection with view cell
-		{
-			VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(node);
-      
-			if (!leaf->mViewCell->Mailed())
-			{
-				//ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
-				leaf->mViewCell->Mail();
-				++ hits;
-			}
-			
-   			//-- fetch the next far child from the stack
-			if (tStack.empty())
-				break;
-      
-			entp = extp;
-			mint = maxt; // NOTE: need this?
-
-			if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
-				break;
-
-			VspBspRayTraversalData &s = tStack.top();
-
-			node = s.mNode;
-			extp = s.mExitPoint;
-			maxt = s.mMaxT;
-
-			tStack.pop();
-		}
-	}
-
-	return hits;
-}
-
-bool VspBspTree::Export(const string filename)
-{
-	Exporter *exporter = Exporter::GetExporter(filename);
-
-	if (exporter) 
-	{
-		//exporter->ExportVspBspTree(*this);
-		return true;
-	}	
-
-	return false;
-}
-
-void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
-{
-	stack<VspBspNode *> nodeStack;
-	nodeStack.push(mRoot);
-
-	ViewCell::NewMail();
-
-	while (!nodeStack.empty()) 
-	{
-		VspBspNode *node = nodeStack.top();
-		nodeStack.pop();
-
-		if (node->IsLeaf()) 
-		{
-			ViewCell *viewCell = dynamic_cast<VspBspLeaf *>(node)->mViewCell;
-
-			if (!viewCell->Mailed()) 
-			{
-				viewCell->Mail();
-				viewCells.push_back(viewCell);
-			}
-		}
-		else 
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-			nodeStack.push(interior->mFront);
-			nodeStack.push(interior->mBack);
-		}
-	}
-}
-
-void VspBspTree::EvaluateViewCellsStats(VspBspViewCellsStatistics &stat) const
-{
-	stat.Reset();
-
-	stack<VspBspNode *> nodeStack;
-	nodeStack.push(mRoot);
-
-	ViewCell::NewMail();
-
-	// exclude root cell
-	mRootCell->Mail();
-
-	while (!nodeStack.empty()) 
-	{
-		VspBspNode *node = nodeStack.top();
-		nodeStack.pop();
-
-		if (node->IsLeaf()) 
-		{
-			++ stat.bspLeaves;
-
-			BspViewCell *viewCell = dynamic_cast<VspBspLeaf *>(node)->mViewCell;
-
-			if (!viewCell->Mailed()) 
-			{
-				viewCell->Mail();
-				
-				++ stat.viewCells;
-				const int pvsSize = viewCell->GetPvs().GetSize();
-
-                stat.pvs += pvsSize;
-
-				if (pvsSize < 1)
-					++ stat.emptyPvs;
-
-				if (pvsSize > stat.maxPvs)
-					stat.maxPvs = pvsSize;
-
-				if (pvsSize < stat.minPvs)
-					stat.minPvs = pvsSize;
-
-				//TODO
-				//if ((int)viewCell->mVspBspLeaves.size() > stat.maxVspBspLeaves)
-				//	stat.maxVspBspLeaves = (int)viewCell->mVspBspLeaves.size();		
-			}
-		}
-		else 
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-			nodeStack.push(interior->mFront);
-			nodeStack.push(interior->mBack);
-		}
-	}
-}
-
-VspBspTreeStatistics &VspBspTree::GetStat()
-{
-	return mStat;
-}
-
-float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
-{
-	float len = 0;
-
-	RayInfoContainer::const_iterator it, it_end = rays.end();
-
-	for (it = rays.begin(); it != it_end; ++ it)
-		len += (*it).SegmentLength();
-
-	return len;
-}
-
-int VspBspTree::SplitRays(const Plane3 &plane,
-					   RayInfoContainer &rays, 
-					   RayInfoContainer &frontRays, 
-					   RayInfoContainer &backRays)
-{
-	int splits = 0;
-
-	while (!rays.empty())
-	{
-		RayInfo bRay = rays.back();
-		VssRay *ray = bRay.mRay;
-
-		rays.pop_back();
-	
-		// get classification and receive new t
-		const int cf = 	bRay.ComputeRayIntersection(plane, ray->mT);
-
-		switch(cf)
-		{
-		case -1:
-			backRays.push_back(bRay);
-			break;
-		case 1:
-			break;
-			frontRays.push_back(bRay);
-		case 0: 
-			//-- split ray
-
-			// if start point behind plane
-			if (plane.Side(bRay.ExtrapOrigin()) <= 0)
-			{
-				backRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
-				frontRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
-			}
-			else
-			{
-				frontRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
-				backRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
-			}
-			break;
-		default:
-			Debug << "Should not come here 4" << endl;
-			break;
-		}
-	}
-
-	return splits;
-}
-
-void VspBspTree::ExtractHalfSpaces(VspBspNode *n, vector<Plane3> &halfSpaces) const
-{
-	VspBspNode *lastNode;
-	do
-	{
-		lastNode = n;
-
-		// want to get planes defining geometry of this node => don't take
-		// split plane of node itself
-		n = n->GetParent();
-		
-		if (n)
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(n);
-			Plane3 halfSpace = *dynamic_cast<VspBspInterior *>(interior)->GetPlane();
-
-            if (interior->mFront != lastNode)
-				halfSpace.ReverseOrientation();
-
-			halfSpaces.push_back(halfSpace);
-		}
-	}
-	while (n);
-}
-
-void VspBspTree::ConstructGeometry(VspBspNode *n, VspBspNodeGeometry &cell) const
-{
-	PolygonContainer polys;
-	ConstructGeometry(n, polys);
-	cell.mPolys = polys;
-}
-
-void VspBspTree::ConstructGeometry(VspBspNode *n, PolygonContainer &cell) const
-{
-	vector<Plane3> halfSpaces;
-	ExtractHalfSpaces(n, halfSpaces);
-
-	PolygonContainer candidatePolys;
-
-	// bounded planes are added to the polygons (reverse polygons
-	// as they have to be outfacing
-	for (int i = 0; i < (int)halfSpaces.size(); ++ i)
-	{
-		Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
-		
-		if (p->Valid())
-		{
-			candidatePolys.push_back(p->CreateReversePolygon());
-			DEL_PTR(p);
-		}
-	}
-
-	// add faces of bounding box (also could be faces of the cell)
-	for (int i = 0; i < 6; ++ i)
-	{
-		VertexContainer vertices;
-	
-		for (int j = 0; j < 4; ++ j)
-			vertices.push_back(mBox.GetFace(i).mVertices[j]);
-
-		candidatePolys.push_back(new Polygon3(vertices));
-	}
-
-	for (int i = 0; i < (int)candidatePolys.size(); ++ i)
-	{
-		// polygon is split by all other planes
-		for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
-		{
-			if (i == j) // polygon and plane are coincident
-				continue;
-
-			VertexContainer splitPts;
-			Polygon3 *frontPoly, *backPoly;
-
-			const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
-			
-			switch (cf)
-			{
-				case Polygon3::SPLIT:
-					frontPoly = new Polygon3();
-					backPoly = new Polygon3();
-
-					candidatePolys[i]->Split(halfSpaces[j], 
-											 *frontPoly, 
-											 *backPoly);
-
-					DEL_PTR(candidatePolys[i]);
-
-					if (frontPoly->Valid())
-						candidatePolys[i] = frontPoly;
-					else
-						DEL_PTR(frontPoly);
-
-					DEL_PTR(backPoly);
-					break;
-				case Polygon3::BACK_SIDE:
-					DEL_PTR(candidatePolys[i]);
-					break;
-				// just take polygon as it is
-				case Polygon3::FRONT_SIDE:
-				case Polygon3::COINCIDENT:
-				default:
-					break;
-			}
-		}
-		
-		if (candidatePolys[i])
-			cell.push_back(candidatePolys[i]);
-	}
-}
-
-void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &cell) const
-{
-	/*
-	TODO
-	vector<VspBspLeaf *> leaves = vc->mBspLeaves;
-
-	vector<VspBspLeaf *>::const_iterator it, it_end = leaves.end();
-
-	for (it = leaves.begin(); it != it_end; ++ it)
-		ConstructGeometry(*it, vcGeom);*/
-}
-
-int VspBspTree::FindNeighbors(VspBspNode *n, vector<VspBspLeaf *> &neighbors, 
-						   const bool onlyUnmailed) const
-{
-	PolygonContainer cell;
-
-	ConstructGeometry(n, cell);
-
-	stack<VspBspNode *> nodeStack;
-	nodeStack.push(mRoot);
-		
-	// planes needed to verify that we found neighbor leaf.
-	vector<Plane3> halfSpaces;
-	ExtractHalfSpaces(n, halfSpaces);
-
-	while (!nodeStack.empty()) 
-	{
-		VspBspNode *node = nodeStack.top();
-		nodeStack.pop();
-
-		if (node->IsLeaf())
-		{
-            if (node != n && (!onlyUnmailed || !node->Mailed()))
-			{
-				// test all planes of current node if candidate really
-				// is neighbour
-				PolygonContainer neighborCandidate;
-				ConstructGeometry(node, neighborCandidate);
-				
-				bool isAdjacent = true;
-				for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
-				{
-					const int cf = 
-						Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);
-
-					if (cf == Polygon3::BACK_SIDE)
-						isAdjacent = false;
-				}
-
-				if (isAdjacent)
-					neighbors.push_back(dynamic_cast<VspBspLeaf *>(node));
-
-				CLEAR_CONTAINER(neighborCandidate);
-			}
-		}
-		else 
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-	
-			const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);
-
-			if (cf == Polygon3::FRONT_SIDE)
-				nodeStack.push(interior->mFront);
-			else
-				if (cf == Polygon3::BACK_SIDE)
-					nodeStack.push(interior->mBack);
-				else 
-				{
-					// random decision
-					nodeStack.push(interior->mBack);
-					nodeStack.push(interior->mFront);
-				}
-		}
-	}
-	
-	CLEAR_CONTAINER(cell);
-	return (int)neighbors.size();
-}
-
-VspBspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
-{
-    stack<VspBspNode *> nodeStack;
-  	nodeStack.push(mRoot);
-	
-	int mask = rand();
-  
-	while (!nodeStack.empty()) 
-	{
-		VspBspNode *node = nodeStack.top();
-		nodeStack.pop();
-	  
-		if (node->IsLeaf()) 
-		{
-			return dynamic_cast<VspBspLeaf *>(node);
-		} 
-		else 
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-			
-			VspBspNode *next;
-	
-			PolygonContainer cell;
-
-			// todo: not very efficient: constructs full cell everytime
-			ConstructGeometry(interior, cell);
-
-			const int cf = Polygon3::ClassifyPlane(cell, halfspace);
-
-			if (cf == Polygon3::BACK_SIDE)
-				next = interior->mFront;
-			else
-				if (cf == Polygon3::FRONT_SIDE)
-					next = interior->mFront;
-			else 
-			{
-				// random decision
-				if (mask & 1)
-					next = interior->mBack;
-				else
-					next = interior->mFront;
-				mask = mask >> 1;
-			}
-
-			nodeStack.push(next);
-		}
-	}
-	
-	return NULL;
-}
-
-VspBspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
-{
-	stack<VspBspNode *> nodeStack;
-	
-	nodeStack.push(mRoot);
-
-	int mask = rand();
-	
-	while (!nodeStack.empty()) 
-	{
-		VspBspNode *node = nodeStack.top();
-		nodeStack.pop();
-		
-		if (node->IsLeaf()) 
-		{
-			if ( (!onlyUnmailed || !node->Mailed()) )
-				return dynamic_cast<VspBspLeaf *>(node);
-		}
-		else 
-		{
-			VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-			// random decision
-			if (mask & 1)
-				nodeStack.push(interior->mBack);
-			else
-				nodeStack.push(interior->mFront);
-
-			mask = mask >> 1;
-		}
-	}
-	
-	return NULL;
-}
-
-int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
-{
-	int pvsSize = 0;
-
-	RayInfoContainer::const_iterator rit, rit_end = rays.end();
-
-	Intersectable::NewMail();
-
-	for (rit = rays.begin(); rit != rays.end(); ++ rit)
-	{
-		VssRay *ray = (*rit).mRay;
-		
-		if (ray->mOriginObject)
-		{
-			if (!ray->mOriginObject->Mailed())
-			{
-				ray->mOriginObject->Mail();
-				++ pvsSize;
-			}
-		}
-		if (ray->mTerminationObject)
-		{
-			if (!ray->mTerminationObject->Mailed())
-			{
-				ray->mTerminationObject->Mail();
-				++ pvsSize;
-			}
-		}
-	}
-
-	return pvsSize;
-}
-
-/*************************************************************
- *            VspBspNodeGeometry Implementation                 *
- *************************************************************/
-
-VspBspNodeGeometry::~VspBspNodeGeometry()
-{
-	CLEAR_CONTAINER(mPolys);
-}
-
-float VspBspNodeGeometry::GetArea() const 
-{ 
-	return Polygon3::GetArea(mPolys);
-}
-
-void VspBspNodeGeometry::SplitGeometry(VspBspNodeGeometry &front,
-									VspBspNodeGeometry &back,
-									const VspBspTree &tree,						 
-									const Plane3 &splitPlane) const
-{	
-	// get cross section of new polygon
-	Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);
-
-	planePoly = SplitPolygon(planePoly, tree);
-
-	//-- plane poly splits all other cell polygons
-	for (int i = 0; i < (int)mPolys.size(); ++ i)
-	{
-		const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
-			
-		// split new polygon with all previous planes
-		switch (cf)
-		{
-			case Polygon3::SPLIT:
-				{
-					Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);
-
-					Polygon3 *frontPoly = new Polygon3();
-					Polygon3 *backPoly = new Polygon3();
-				
-					poly->Split(splitPlane, *frontPoly, *backPoly);
-
-					DEL_PTR(poly);
-
-					if (frontPoly->Valid())
-						front.mPolys.push_back(frontPoly);
-					if (backPoly->Valid())
-						back.mPolys.push_back(backPoly);
-				}
-				
-				break;
-			case Polygon3::BACK_SIDE:
-				back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));			
-				break;
-			case Polygon3::FRONT_SIDE:
-				front.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));	
-				break;
-			case Polygon3::COINCIDENT:
-				//front.mPolys.push_back(CreateReversePolygon(mPolys[i]));
-				back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));
-				break;
-			default:
-				break;
-		}
-	}
-
-	//-- finally add the new polygon to the child cells
-	if (planePoly)
-	{
-		// add polygon with normal pointing into positive half space to back cell
-		back.mPolys.push_back(planePoly);
-		// add polygon with reverse orientation to front cell
-		front.mPolys.push_back(planePoly->CreateReversePolygon());
-	}
-
-	//Debug << "returning new geometry " << mPolys.size() << " f: " << front.mPolys.size() << " b: " << back.mPolys.size() << endl;
-	//Debug << "old area " << GetArea() << " f: " << front.GetArea() << " b: " << back.GetArea() << endl;
-}
-
-Polygon3 *VspBspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
-										const VspBspTree &tree) const
-{
-	// polygon is split by all other planes
-	for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
-	{
-		Plane3 plane = mPolys[i]->GetSupportingPlane();
-
-		const int cf = 
-			planePoly->ClassifyPlane(plane, 0.00001f);
-			
-		// split new polygon with all previous planes
-		switch (cf)
-		{
-			case Polygon3::SPLIT:
-				{
-					Polygon3 *frontPoly = new Polygon3();
-					Polygon3 *backPoly = new Polygon3();
-
-					planePoly->Split(plane, *frontPoly, *backPoly);
-					DEL_PTR(planePoly);
-
-					if (backPoly->Valid())
-						planePoly = backPoly;
-					else
-						DEL_PTR(backPoly);
-				}
-				break;
-			case Polygon3::FRONT_SIDE:
-				DEL_PTR(planePoly);
-                break;
-			// polygon is taken as it is
-			case Polygon3::BACK_SIDE:
-			case Polygon3::COINCIDENT:
-			default:
-				break;
-		}
-	}
-
-	return planePoly;
-}
-
-void VspBspViewCellsStatistics::Print(ostream &app) const
-{
-	app << "===== VspBspViewCells statistics ===============\n";
-
-	app << setprecision(4);
-
-	//app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
-
-	app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
-
-	app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
-
-	app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
-
-	app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
-
-	app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
-
-	app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
-
-	app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgVspBspLeaves() << endl;
-
-	app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxVspBspLeaves << endl;
-	
-	app << "===== END OF VspBspViewCells statistics ==========\n";
-}
-
-BspViewCell *VspBspTree::GetRootCell() const
-{
-	return mRootCell;
-}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 428)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 448)
@@ -15,4 +15,5 @@
 class VssRay;
 
+
 // -------------------------------------------------------------------
 // CRay class.  A ray is defined by a location and a direction.
@@ -48,5 +49,6 @@
 		
     bool operator<(
-									 const Intersection &b) const {
+				   const Intersection &b) const {
+
       return 
 				mT
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 446)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp	(revision 448)
@@ -1269,5 +1269,5 @@
 	  "100000");
 
-  RegisterOption("BspTree.Construction.sideTolerance",
+  RegisterOption("BspTree.Construction.epsilon",
 	  optFloat,
 	  "-bsp_construction_side_tolerance=",
@@ -1732,7 +1732,17 @@
 	  "false");
 
-   RegisterOption("VspBspTree.Factor.leastRaySplits", optFloat, "-vsp_bsp_factor_least_ray_splits=", "1.0");
-   RegisterOption("VspBspTree.Factor.balancedRays", optFloat, "-vsp_bsp_factor_balanced_rays=", "1.0");
-   RegisterOption("VspBspTree.Factor.pvs", optFloat, "-vsp_bsp_factor_pvs=", "1.0");
+  RegisterOption("VspBspTree.Construction.samples",
+	  optInt,
+	  "-bsp_construction_samples=",
+	  "100000");
+
+  RegisterOption("VspBspTree.Construction.epsilon",
+	  optFloat,
+	  "-vsp_bsp_construction_side_tolerance=",
+	  "0.002");
+
+  RegisterOption("VspBspTree.Factor.leastRaySplits", optFloat, "-vsp_bsp_factor_least_ray_splits=", "1.0");
+  RegisterOption("VspBspTree.Factor.balancedRays", optFloat, "-vsp_bsp_factor_balanced_rays=", "1.0");
+  RegisterOption("VspBspTree.Factor.pvs", optFloat, "-vsp_bsp_factor_pvs=", "1.0");
 
    //////////////////////////////////////////////////////////////////////////////////
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h	(revision 446)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h	(revision 448)
@@ -20,4 +20,5 @@
 class Polygon3;
 class VssTree;
+class VspBspTree;
 class RssTree;
 
@@ -44,5 +45,4 @@
   virtual bool ExportScene(SceneGraphNode *node) = 0;
 
-
   virtual bool
   ExportBox(const AxisAlignedBox3 &box) = 0;
@@ -56,9 +56,4 @@
   virtual bool
   ExportVssTree2(const VssTree &tree,
-				 const Vector3 direction
-				 ) = 0;
-
-  virtual bool
-  ExportRssTree2(const RssTree &tree,
 				 const Vector3 direction
 				 ) = 0;
@@ -98,7 +93,14 @@
   virtual void 
   ExportBspSplits(const BspTree &tree, const bool exportDepth = false) = 0;
+
   virtual void
   ExportLeavesGeometry(const BspTree &tree, const vector<BspLeaf *> &leaves) = 0;
-	
+  
+  virtual void 
+  ExportBspSplits(const VspBspTree &tree, const bool exportDepth = false) = 0;
+
+  virtual void
+  ExportBspViewCellPartition(const VspBspTree &tree, const int maxPvs = 0) = 0;
+
   virtual void 
   ExportPolygons(const PolygonContainer &polys) = 0;
@@ -112,4 +114,8 @@
   virtual void
   ExportGeometry(const ObjectContainer &objects) = 0;
+
+  virtual bool
+  ExportRssTree2(const RssTree &tree,
+				 const Vector3 direction) = 0;
 
   void SetExportRayDensity(const bool d) { mExportRayDensity = d; }
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp	(revision 440)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp	(revision 448)
@@ -3,4 +3,6 @@
 #include "ViewCellBsp.h"
 #include "ViewCell.h"
+#include "VspBspTree.h"
+
 
 void SimulationStatistics::Print(ostream &app) const
@@ -53,4 +55,6 @@
 /*     class BspRenderSimulator implementation          */
 /********************************************************/
+
+
 BspRenderSimulator::BspRenderSimulator(BspTree *bspTree):
 mBspTree(bspTree)
@@ -139,5 +143,5 @@
 
 Real BspRenderSimulator::RenderPvs(ViewCell &viewCell, 
-									   float objRenderTime) const
+								   float objRenderTime) const
 {
 	return viewCell.GetPvs().GetSize() * objRenderTime;
@@ -234,2 +238,98 @@
 	return leaf->mKdPvs.GetSize() * objRenderTime;
 }
+
+/********************************************************/
+/*     class BspRenderSimulator implementation          */
+/********************************************************/
+VspBspRenderSimulator::VspBspRenderSimulator(VspBspTree *vspBspTree):
+mVspBspTree(vspBspTree)
+{
+}
+
+VspBspRenderSimulator::VspBspRenderSimulator(float objRenderCost, 
+											 float vcOverhead, 
+											 float moveSpeed,
+											 VspBspTree *vspBspTree):
+RenderSimulator(objRenderCost, vcOverhead, moveSpeed), 
+mVspBspTree(vspBspTree)
+{
+}
+
+SimulationStatistics VspBspRenderSimulator::SimulateRendering()
+{
+	SimulationStatistics simStat;
+
+	simStat.Reset();
+	simStat.Start();
+
+	Real renderTime = 0;
+	
+	// overhead for loading the PVS of the view cells
+	float loadPvsOverhead = 0;
+	// probability that view point lies in a view cell
+	float pInVcTotal = 0;
+
+	// total probability that a view cell border is crossed
+	const float pCrossVcTotal = mVspBspTree->GetBoundingBox().SurfaceArea();
+
+	// collect all view cells
+	ViewCellContainer viewCells;
+	mVspBspTree->CollectViewCells(viewCells);
+
+	ViewCellContainer::const_iterator it, it_end = viewCells.end();
+
+	// surface area substitute for probability
+	PolygonContainer geom;
+float overallarea = 0;
+	for (it = viewCells.begin(); it != it_end; ++ it)
+	{
+		// compute view cell area
+		mVspBspTree->ConstructGeometry(dynamic_cast<BspViewCell *>(*it), geom);
+
+		const float area = Polygon3::GetArea(geom);
+		if (area < 0.0001)
+			Debug << "warning, area: " << area << endl;
+		CLEAR_CONTAINER(geom);
+
+		// area substitute for view point probability
+		float pInVc = area;
+				
+		// compute render time of PVS times probability that view point is in view cell
+		float vcCost = pInVc * RenderPvs(*(*it), mObjRenderCost);
+		//Debug << "p: " << pInVc << " rendercost: " << RenderPvs(*(*it), mObjRenderCost) << endl;
+		renderTime += vcCost;
+
+		if (vcCost > simStat.maxCost)
+			simStat.maxCost = vcCost;
+		else if (vcCost < simStat.minCost)
+			simStat.minCost = vcCost;
+
+		// probability that a view cell border is crossed
+		float pCrossVc = area * mMoveSpeed;
+
+		// crossing the border of a view cell is also depending on the move speed
+		loadPvsOverhead += pCrossVc * mVcOverhead;
+overallarea += area;
+		pInVcTotal += pInVc;
+	}
+	Debug << "overall area: " << overallarea << endl;
+	
+	renderTime /= pInVcTotal;
+	loadPvsOverhead /= pCrossVcTotal;
+
+	simStat.avgRtWithoutOverhead = renderTime;
+	simStat.avgRenderTime = renderTime + loadPvsOverhead;
+	
+	simStat.maxCost /= pInVcTotal;
+	simStat.minCost /= pInVcTotal;
+
+	simStat.Stop();
+
+	return simStat;
+}
+
+Real VspBspRenderSimulator::RenderPvs(ViewCell &viewCell, 
+									  float objRenderTime) const
+{
+	return viewCell.GetPvs().GetSize() * objRenderTime;
+}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 445)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 448)
@@ -98,7 +98,7 @@
 }
 
-Plane3 *BspInterior::GetPlane()
-{
-	return &mPlane;
+Plane3 BspInterior::GetPlane() const
+{
+	return mPlane;
 }
 
@@ -117,81 +117,13 @@
 }
 
-int BspInterior::SplitPolygons(PolygonContainer &polys, 
-							   PolygonContainer &frontPolys, 
-							   PolygonContainer &backPolys, 
-							   PolygonContainer &coincident)
-{
-	int splits = 0;
-
-#ifdef _Debug
-	Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
-#endif
-	while (!polys.empty())
-	{
-		Polygon3 *poly = polys.back();
-		polys.pop_back();
-
-		//Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";
-
-		// classify polygon
-		const int cf = poly->ClassifyPlane(mPlane);
-
-		switch (cf)
-		{
-			case Polygon3::COINCIDENT:
-				coincident.push_back(poly);
-				break;			
-			case Polygon3::FRONT_SIDE:	
-				frontPolys.push_back(poly);
-				break;
-			case Polygon3::BACK_SIDE:
-				backPolys.push_back(poly);
-				break;
-			case Polygon3::SPLIT:
-				{
-					Polygon3 *front_piece = new Polygon3(poly->mParent);
-					Polygon3 *back_piece = new Polygon3(poly->mParent);
-
-					//-- split polygon into front and back part
-					poly->Split(mPlane, *front_piece, *back_piece);
-					
-					++ splits; // increase number of splits
-
-					//-- inherit rays from parent polygon for blocked ray criterium
-					poly->InheritRays(*front_piece, *back_piece);
-					//Debug << "p: " << poly->mPiercingRays.size() << " f: " << front_piece->mPiercingRays.size() << " b: " << back_piece->mPiercingRays.size() << endl;
-				
-					// check if polygons still valid 
-					if (front_piece->Valid())
-						frontPolys.push_back(front_piece);
-					else
-						DEL_PTR(front_piece);
-				
-					if (back_piece->Valid())
-						backPolys.push_back(back_piece);
-					else				
-						DEL_PTR(back_piece);
-				
-#ifdef _DEBUG
-					Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
-#endif
-					DEL_PTR(poly);
-				}
-				break;
-			default:
-                Debug << "SHOULD NEVER COME HERE\n";
-				break;
-		}
-	}
-
-	return splits;
-}
-
 /****************************************************************/
 /*                  class BspLeaf implementation                */
 /****************************************************************/
+
+
 BspLeaf::BspLeaf(): mViewCell(NULL)
 {
 }
+
 
 BspLeaf::BspLeaf(BspViewCell *viewCell): 
@@ -199,4 +131,5 @@
 {
 }
+
 
 BspLeaf::BspLeaf(BspInterior *parent): 
@@ -204,4 +137,6 @@
 {}
 
+
+
 BspLeaf::BspLeaf(BspInterior *parent, BspViewCell *viewCell): 
 BspNode(parent), mViewCell(viewCell)
@@ -222,36 +157,4 @@
 { 
 	return true; 
-}
-
-void BspLeaf::AddToPvs(const BoundedRayContainer &rays, 
-					   int &sampleContributions,
-					   int &contributingSamples)
-{
-	sampleContributions = 0;
-	contributingSamples = 0;
-
-    BoundedRayContainer::const_iterator it, it_end = rays.end();
-
-	// add contributions from samples to the PVS
-	for (it = rays.begin(); it != it_end; ++ it)
-	{
-		int contribution = 0;
-		Ray *ray = (*it)->mRay;
-			
-		if (!ray->intersections.empty())
-			contribution += mViewCell->GetPvs().AddSample(ray->intersections[0].mObject);
-		
-		if (ray->sourceObject.mObject)
-			contribution += mViewCell->GetPvs().AddSample(ray->sourceObject.mObject);
-
-		if (contribution)
-		{
-			sampleContributions += contribution;
-			++ contributingSamples;
-		}
-
-		//if (ray->mFlags & Ray::STORE_BSP_INTERSECTIONS)
-			ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
-	}
 }
 
@@ -308,8 +211,6 @@
 	environment->GetIntValue("BspTree.maxTests", mMaxTests);
 
-	environment->GetFloatValue("BspTree.Construction.sideTolerance", Vector3::sDistTolerance);
-	Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
-
-
+	environment->GetFloatValue("BspTree.Construction.epsilon", mEpsilon);
+	
     Debug << "BSP max depth: " << mTermMaxDepth << endl;
 	Debug << "BSP min PVS: " << mTermMinPvs << endl;
@@ -347,8 +248,85 @@
 }
 
+
 const BspTreeStatistics &BspTree::GetStatistics() const 
 {
 	return mStat;
 }
+
+
+int BspTree::SplitPolygons(const Plane3 &plane,
+						   PolygonContainer &polys, 
+						   PolygonContainer &frontPolys, 
+						   PolygonContainer &backPolys, 
+						   PolygonContainer &coincident) const
+{
+	int splits = 0;
+
+#ifdef _Debug
+	Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
+#endif
+	while (!polys.empty())
+	{
+		Polygon3 *poly = polys.back();
+		polys.pop_back();
+
+		//Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";
+
+		// classify polygon
+		const int cf = poly->ClassifyPlane(plane, mEpsilon);
+
+		switch (cf)
+		{
+			case Polygon3::COINCIDENT:
+				coincident.push_back(poly);
+				break;			
+			case Polygon3::FRONT_SIDE:	
+				frontPolys.push_back(poly);
+				break;
+			case Polygon3::BACK_SIDE:
+				backPolys.push_back(poly);
+				break;
+			case Polygon3::SPLIT:
+				{
+					Polygon3 *front_piece = new Polygon3(poly->mParent);
+					Polygon3 *back_piece = new Polygon3(poly->mParent);
+
+					//-- split polygon into front and back part
+					poly->Split(plane, 
+								*front_piece, 
+								*back_piece,
+								mEpsilon);
+					
+					++ splits; // increase number of splits
+
+					//-- inherit rays from parent polygon for blocked ray criterium
+					poly->InheritRays(*front_piece, *back_piece);
+				
+					// check if polygons still valid 
+					if (front_piece->Valid(mEpsilon))
+						frontPolys.push_back(front_piece);
+					else
+						DEL_PTR(front_piece);
+				
+					if (back_piece->Valid(mEpsilon))
+						backPolys.push_back(back_piece);
+					else				
+						DEL_PTR(back_piece);
+				
+#ifdef _DEBUG
+					Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
+#endif
+					DEL_PTR(poly);
+				}
+				break;
+			default:
+                Debug << "SHOULD NEVER COME HERE\n";
+				break;
+		}
+	}
+
+	return splits;
+}
+
 
 void BspTreeStatistics::Print(ostream &app) const
@@ -458,8 +436,9 @@
 		
 				// split viewcell polygons with respect to split plane
-				splits += interior->SplitPolygons(*tData.mPolygons,
-												  *frontPolys, 
-												  *backPolys,
-												  coincident);
+				splits += SplitPolygons(interior->GetPlane(),
+										*tData.mPolygons, 
+										*frontPolys, 
+										*backPolys, 
+										coincident);
 				
 				// extract view cells associated with the split polygons
@@ -505,4 +484,5 @@
 			// cleanup
 			DEL_PTR(tData.mPolygons);
+			DEL_PTR(tData.mRays);
 		}
 		else
@@ -514,5 +494,7 @@
 }
 
-int BspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
+int BspTree::AddMeshToPolygons(Mesh *mesh, 
+							   PolygonContainer &polys, 
+							   MeshInstance *parent)
 {
 	FaceContainer::const_iterator fi;
@@ -523,5 +505,5 @@
 		Polygon3 *poly = new Polygon3((*fi), mesh);
 		
-		if (poly->Valid())
+		if (poly->Valid(mEpsilon))
 		{
 			poly->mParent = parent; // set parent intersectable
@@ -684,5 +666,5 @@
 
 		float minT, maxT;
-		if (BoundRay(*ray, minT, maxT))
+		if (mBox.GetRaySegment(*ray, minT, maxT))
 			rays->push_back(new BoundedRay(ray, minT, maxT));
 	}
@@ -719,5 +701,5 @@
 
 		float minT, maxT;
-		if (BoundRay(*ray, minT, maxT))
+		if (mBox.GetRaySegment(*ray, minT, maxT))
 			rays->push_back(new BoundedRay(ray, minT, maxT));
 	}
@@ -801,5 +783,5 @@
 		{
 			int conSamp = 0, sampCon = 0;
-			leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
+			AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
 			
 			mStat.contributingSamples += conSamp;
@@ -818,5 +800,6 @@
 		DEL_PTR(tData.mPolygons);
 		DEL_PTR(tData.mRays);
-			
+		DEL_PTR(tData.mGeometry);
+
 		return leaf;
 	}
@@ -860,4 +843,5 @@
 	// cleanup
 	DEL_PTR(tData.mNode);
+
 	DEL_PTR(tData.mPolygons);
 	DEL_PTR(tData.mRays);
@@ -916,6 +900,6 @@
 	
 
-	Debug << "number of rays: " << tData.mRays->size() << endl;
-	Debug << "number of polys: " << tData.mPolygons->size() << endl;
+	Debug << "number of rays: " << (int)tData.mRays->size() << endl;
+	Debug << "number of polys: " << (int)tData.mPolygons->size() << endl;
 
 	startTime = GetTime();
@@ -928,8 +912,9 @@
 	startTime = GetTime();
 	// subdivide polygons with plane
-	mStat.splits += interior->SplitPolygons(*tData.mPolygons, 
-		                                    *frontData.mPolygons, 
-											*backData.mPolygons, 
-											coincident);
+	mStat.splits += SplitPolygons(interior->GetPlane(),
+								  *tData.mPolygons, 
+		                          *frontData.mPolygons, 
+								  *backData.mPolygons, 
+								  coincident);
 
 	Debug << "time used for polygon splitting: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
@@ -944,6 +929,7 @@
 		tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
 									   *backData.mGeometry, 
-									   *this, 
-									   interior->mPlane);
+									   interior->mPlane,
+									   mBox,
+									   mEpsilon);
 	
 		
@@ -974,6 +960,6 @@
 	interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
 	
-	frontData.mNode = interior->mFront;
-	backData.mNode = interior->mBack;
+	frontData.mNode = interior->GetFront();
+	backData.mNode = interior->GetBack();
 	
 	//DEL_PTR(leaf);
@@ -1269,5 +1255,6 @@
 
 			BoundedRay *bRay = (*rays)[idx[j]];
-			pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : bRay->mRay->Extrap(bRay->mMaxT);
+			pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : 
+								bRay->mRay->Extrap(bRay->mMaxT);
 		}	
 			
@@ -1347,5 +1334,6 @@
 		Polygon3 *poly = polys[testIdx];
 
-        const int classification = poly->ClassifyPlane(candidatePlane);
+        const int classification = 
+			poly->ClassifyPlane(candidatePlane, mEpsilon);
 
 		if (mSplitPlaneStrategy & BALANCED_POLYS)
@@ -1437,26 +1425,4 @@
 }
 
-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;
-}
 
 inline void BspTree::GenerateUniqueIdsForPvs()
@@ -1499,5 +1465,9 @@
 			BspNodeGeometry backCell;
 
-			cell.SplitGeometry(frontCell, backCell, *this, candidatePlane);
+			cell.SplitGeometry(frontCell, 
+							   backCell, 
+							   candidatePlane,
+							   mBox,
+							   mEpsilon);
 		
 			pFront = frontCell.GetArea();
@@ -1803,5 +1773,5 @@
 	float maxt, mint;
 
-	if (!BoundRay(ray, mint, maxt))
+	if (!mBox.GetRaySegment(ray, mint, maxt))
 		return 0;
 
@@ -1818,12 +1788,9 @@
 		if (!node->IsLeaf()) 
 		{
-			BspInterior *in = (BspInterior *) node;
+			BspInterior *in = dynamic_cast<BspInterior *>(node);
 			
-			Plane3 *splitPlane = in->GetPlane();
-
-			int entSide = splitPlane->Side(entp);
-			int extSide = splitPlane->Side(extp);
-
-			Vector3 intersection;
+			Plane3 splitPlane = in->GetPlane();
+			const int entSide = splitPlane.Side(entp);
+			const int extSide = splitPlane.Side(extp);
 
 			if (entSide < 0)
@@ -1858,5 +1825,5 @@
 			// find intersection of ray segment with plane
 			float t;
-			extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &t);
+			extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
 			maxt *= t;
 			
@@ -1934,6 +1901,6 @@
 			BspInterior *interior = dynamic_cast<BspInterior *>(node);
 
-			nodeStack.push(interior->mFront);
-			nodeStack.push(interior->mBack);
+			nodeStack.push(interior->GetFront());
+			nodeStack.push(interior->GetBack());
 		}
 	}
@@ -1989,6 +1956,6 @@
 			BspInterior *interior = dynamic_cast<BspInterior *>(node);
 
-			nodeStack.push(interior->mFront);
-			nodeStack.push(interior->mBack);
+			nodeStack.push(interior->GetFront());
+			nodeStack.push(interior->GetBack());
 		}
 	}
@@ -2101,7 +2068,7 @@
 		{
 			BspInterior *interior = dynamic_cast<BspInterior *>(n);
-			Plane3 halfSpace = *dynamic_cast<BspInterior *>(interior)->GetPlane();
-
-            if (interior->mFront != lastNode)
+			Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
+
+            if (interior->GetFront() != lastNode)
 				halfSpace.ReverseOrientation();
 
@@ -2135,5 +2102,5 @@
 	ExtractHalfSpaces(n, halfSpaces);
 
-	PolygonContainer candidatePolys;
+	PolygonContainer candidates;
 
 	// bounded planes are added to the polygons (reverse polygons
@@ -2143,7 +2110,7 @@
 		Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
 		
-		if (p->Valid())
-		{
-			candidatePolys.push_back(p->CreateReversePolygon());
+		if (p->Valid(mEpsilon))
+		{
+			candidates.push_back(p->CreateReversePolygon());
 			DEL_PTR(p);
 		}
@@ -2158,11 +2125,11 @@
 			vertices.push_back(mBox.GetFace(i).mVertices[j]);
 
-		candidatePolys.push_back(new Polygon3(vertices));
-	}
-
-	for (int i = 0; i < (int)candidatePolys.size(); ++ i)
+		candidates.push_back(new Polygon3(vertices));
+	}
+
+	for (int i = 0; i < (int)candidates.size(); ++ i)
 	{
 		// polygon is split by all other planes
-		for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
+		for (int j = 0; (j < (int)halfSpaces.size()) && candidates[i]; ++ j)
 		{
 			if (i == j) // polygon and plane are coincident
@@ -2172,5 +2139,6 @@
 			Polygon3 *frontPoly, *backPoly;
 
-			const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
+			const int cf = candidates[i]->
+				ClassifyPlane(halfSpaces[j], mEpsilon);
 			
 			switch (cf)
@@ -2180,12 +2148,13 @@
 					backPoly = new Polygon3();
 
-					candidatePolys[i]->Split(halfSpaces[j], 
-											 *frontPoly, 
-											 *backPoly);
-
-					DEL_PTR(candidatePolys[i]);
-
-					if (frontPoly->Valid())
-						candidatePolys[i] = frontPoly;
+					candidates[i]->Split(halfSpaces[j], 
+										 *frontPoly, 
+										 *backPoly,
+										 mEpsilon);
+
+					DEL_PTR(candidates[i]);
+
+					if (frontPoly->Valid(mEpsilon))
+						candidates[i] = frontPoly;
 					else
 						DEL_PTR(frontPoly);
@@ -2194,5 +2163,5 @@
 					break;
 				case Polygon3::BACK_SIDE:
-					DEL_PTR(candidatePolys[i]);
+					DEL_PTR(candidates[i]);
 					break;
 				// just take polygon as it is
@@ -2204,6 +2173,6 @@
 		}
 		
-		if (candidatePolys[i])
-			cell.push_back(candidatePolys[i]);
+		if (candidates[i])
+			cell.push_back(candidates[i]);
 	}
 }
@@ -2242,5 +2211,7 @@
 				{
 					const int cf = 
-						Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);
+						Polygon3::ClassifyPlane(neighborCandidate, 
+												halfSpaces[i],
+												mEpsilon);
 
 					if (cf == Polygon3::BACK_SIDE)
@@ -2258,16 +2229,18 @@
 			BspInterior *interior = dynamic_cast<BspInterior *>(node);
 	
-			const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);
+			const int cf = Polygon3::ClassifyPlane(cell, 
+												   interior->mPlane,
+												   mEpsilon);
 
 			if (cf == Polygon3::FRONT_SIDE)
-				nodeStack.push(interior->mFront);
+				nodeStack.push(interior->GetFront());
 			else
 				if (cf == Polygon3::BACK_SIDE)
-					nodeStack.push(interior->mBack);
+					nodeStack.push(interior->GetBack());
 				else 
 				{
 					// random decision
-					nodeStack.push(interior->mBack);
-					nodeStack.push(interior->mFront);
+					nodeStack.push(interior->GetBack());
+					nodeStack.push(interior->GetFront());
 				}
 		}
@@ -2305,18 +2278,20 @@
 			ConstructGeometry(interior, cell);
 
-			const int cf = Polygon3::ClassifyPlane(cell, halfspace);
+			const int cf = Polygon3::ClassifyPlane(cell, 
+												   halfspace, 
+												   mEpsilon);
 
 			if (cf == Polygon3::BACK_SIDE)
-				next = interior->mFront;
+				next = interior->GetFront();
 			else
 				if (cf == Polygon3::FRONT_SIDE)
-					next = interior->mFront;
+					next = interior->GetFront();
 			else 
 			{
 				// random decision
 				if (mask & 1)
-					next = interior->mBack;
+					next = interior->GetBack();
 				else
-					next = interior->mFront;
+					next = interior->GetFront();
 				mask = mask >> 1;
 			}
@@ -2353,7 +2328,7 @@
 			// random decision
 			if (mask & 1)
-				nodeStack.push(interior->mBack);
+				nodeStack.push(interior->GetBack());
 			else
-				nodeStack.push(interior->mFront);
+				nodeStack.push(interior->GetFront());
 
 			mask = mask >> 1;
@@ -2362,4 +2337,39 @@
 	
 	return NULL;
+}
+
+void BspTree::AddToPvs(BspLeaf *leaf,
+					   const BoundedRayContainer &rays, 
+					   int &sampleContributions,
+					   int &contributingSamples)
+{
+	sampleContributions = 0;
+	contributingSamples = 0;
+
+    BoundedRayContainer::const_iterator it, it_end = rays.end();
+
+	ViewCell *vc = leaf->GetViewCell();
+
+	// add contributions from samples to the PVS
+	for (it = rays.begin(); it != it_end; ++ it)
+	{
+		int contribution = 0;
+		Ray *ray = (*it)->mRay;
+			
+		if (!ray->intersections.empty())
+			contribution += vc->GetPvs().AddSample(ray->intersections[0].mObject);
+		
+		if (ray->sourceObject.mObject)
+			contribution += vc->GetPvs().AddSample(ray->sourceObject.mObject);
+
+		if (contribution)
+		{
+			sampleContributions += contribution;
+			++ contributingSamples;
+		}
+
+		//if (ray->mFlags & Ray::STORE_BSP_INTERSECTIONS)
+		//	ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
+	}
 }
 
@@ -2397,7 +2407,39 @@
 }
 
-/*************************************************************
- *            BspNodeGeometry Implementation                 *
- *************************************************************/
+float BspTree::GetEpsilon() const
+{
+	return mEpsilon;
+}
+
+void BspViewCellsStatistics::Print(ostream &app) const
+{
+	app << "===== BspViewCells statistics ===============\n";
+
+	app << setprecision(4);
+
+	//app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
+
+	app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
+
+	app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
+
+	app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
+
+	app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
+
+	app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
+
+	app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
+
+	app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgBspLeaves() << endl;
+
+	app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxBspLeaves << endl;
+	
+	app << "===== END OF BspViewCells statistics ==========\n";
+}
+
+/*************************************************************/
+/*            BspNodeGeometry Implementation                 */
+/*************************************************************/
 
 BspNodeGeometry::~BspNodeGeometry()
@@ -2413,18 +2455,21 @@
 void BspNodeGeometry::SplitGeometry(BspNodeGeometry &front,
 									BspNodeGeometry &back,
-									const BspTree &tree,						 
-									const Plane3 &splitPlane) const
+									const Plane3 &splitPlane,
+									const AxisAlignedBox3 &box,
+									const float epsilon) const
 {	
 	// get cross section of new polygon
-	Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);
-
-	planePoly = SplitPolygon(planePoly, tree);
-
-	//-- plane poly splits all other cell polygons
+	Polygon3 *planePoly = box.CrossSection(splitPlane);
+
+	// split polygon with all other polygons
+	planePoly = SplitPolygon(planePoly, epsilon);
+
+	//-- new polygon splits all other polygons
 	for (int i = 0; i < (int)mPolys.size(); ++ i)
 	{
-		const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
+		/// don't use epsilon here to get exact split planes
+		const int cf = 
+			mPolys[i]->ClassifyPlane(splitPlane, Limits::Small);
 			
-		// split new polygon with all previous planes
 		switch (cf)
 		{
@@ -2436,14 +2481,17 @@
 					Polygon3 *backPoly = new Polygon3();
 				
-					poly->Split(splitPlane, *frontPoly, *backPoly);
+					poly->Split(splitPlane, 
+								*frontPoly, 
+								*backPoly,
+								epsilon);
 
 					DEL_PTR(poly);
 
-					if (frontPoly->Valid())
+					if (frontPoly->Valid(epsilon))
 						front.mPolys.push_back(frontPoly);
 					else
 						DEL_PTR(frontPoly);
 
-					if (backPoly->Valid())
+					if (backPoly->Valid(epsilon))
 						back.mPolys.push_back(backPoly);
 					else
@@ -2467,5 +2515,5 @@
 	}
 
-	//-- finally add the new polygon to the child cells
+	//-- finally add the new polygon to the child node geometries
 	if (planePoly)
 	{
@@ -2481,6 +2529,9 @@
 
 Polygon3 *BspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
-										const BspTree &tree) const
-{
+										const float epsilon) const
+{
+	if (!planePoly->Valid(epsilon))
+		DEL_PTR(planePoly);
+
 	// polygon is split by all other planes
 	for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
@@ -2488,6 +2539,7 @@
 		Plane3 plane = mPolys[i]->GetSupportingPlane();
 
+		/// don't use epsilon here to get exact split planes
 		const int cf = 
-			planePoly->ClassifyPlane(plane, 0.00001f);
+			planePoly->ClassifyPlane(plane, Limits::Small);
 			
 		// split new polygon with all previous planes
@@ -2499,5 +2551,8 @@
 					Polygon3 *backPoly = new Polygon3();
 
-					planePoly->Split(plane, *frontPoly, *backPoly);
+					planePoly->Split(plane, 
+									 *frontPoly, 
+									 *backPoly,
+									 epsilon);
 					
 					// don't need anymore
@@ -2506,5 +2561,5 @@
 
 					// back polygon is belonging to geometry
-					if (backPoly->Valid())
+					if (backPoly->Valid(epsilon))
 						planePoly = backPoly;
 					else
@@ -2525,29 +2580,2 @@
 	return planePoly;
 }
-
-void BspViewCellsStatistics::Print(ostream &app) const
-{
-	app << "===== BspViewCells statistics ===============\n";
-
-	app << setprecision(4);
-
-	//app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
-
-	app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
-
-	app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
-
-	app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
-
-	app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
-
-	app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
-
-	app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
-
-	app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgBspLeaves() << endl;
-
-	app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxBspLeaves << endl;
-	
-	app << "===== END OF BspViewCells statistics ==========\n";
-}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp	(revision 448)
@@ -2,7 +2,4 @@
 #include "Vector3.h"
 #include "Halton.h"
-
-float Vector3::sDistTolerance = 0.002f;
-float Vector3::sDistToleranceSqrt = 0.000004f;
 
 // Given min a vector to minimize and a candidate vector, replace
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Preprocessor.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Preprocessor.cpp	(revision 445)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Preprocessor.cpp	(revision 448)
@@ -178,14 +178,18 @@
 		mViewCellsManager = new KdViewCellsManager(mKdTree);
 	}
-	if (strcmp(viewCellsStr, "bspTree") == 0)
+	else if (strcmp(viewCellsStr, "bspTree") == 0)
 	{
 		mBspTree = new BspTree();
+
+		Debug << "view cell type: Bsp" << endl;
 
 		environment->GetIntValue("BspTree.Construction.samples", constructionSamples);
 		mViewCellsManager = new BspViewCellsManager(mBspTree, constructionSamples);
 	}
-	if (strcmp(viewCellsStr, "vspBspTree") == 0)
+	else if (strcmp(viewCellsStr, "vspBspTree") == 0)
 	{
 		mVspBspTree = new VspBspTree();
+
+		Debug << "view cell type: VspBsp" << endl;
 
 		environment->GetIntValue("VspBspTree.Construction.samples", constructionSamples);
@@ -204,4 +208,5 @@
 		mBspTree = new BspTree();
 
+		Debug << "view cell type: Bsp" << endl;
 		environment->GetIntValue("BspTree.Construction.samples", constructionSamples);
 		mViewCellsManager = new BspViewCellsManager(mBspTree, constructionSamples);
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Plane3.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Plane3.h	(revision 437)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Plane3.h	(revision 448)
@@ -51,8 +51,7 @@
   */
   Vector3 FindIntersection(const Vector3 &a,
-			   const Vector3 &b,
-			   float *t = NULL,
-			   bool *coplanar = NULL
-			   ) const 
+						   const Vector3 &b,
+						   float *t = NULL,
+						   bool *coplanar = NULL) const 
   {
     const Vector3 v = b - a; // line from A to B
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 444)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 448)
@@ -225,11 +225,13 @@
 mType(LOCAL_RAY)
 {
-	const float dist = Distance(vssRay.mTermination, vssRay.mOrigin);
-
-	if (dist)
-		dir = vssRay.GetDir() / dist;
+	float len = vssRay.Length();
+
+	if (!len)
+		len = Limits::Small;
+
+	dir = vssRay.GetDir() / len;
 
 	if (vssRay.mTerminationObject)
-		intersections.push_back(Intersection(dist, vssRay.mTerminationObject, 0));
+		intersections.push_back(Intersection(len, vssRay.mTerminationObject, 0));
 
 	Precompute();
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h	(revision 446)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h	(revision 448)
@@ -19,4 +19,6 @@
 class ViewCell;
 class BspTree;
+class VspBspTree;
+class BspNode;
 
 class X3dExporter : public Exporter
@@ -92,4 +94,12 @@
   ExportGeometry(const ObjectContainer &objects);
 
+  bool
+  ExportRays(const RayContainer &rays,
+			 const float length=1000,
+			 const RgbColor &color = RgbColor(1,1,1));
+  bool
+  ExportRays(const VssRayContainer &rays,
+			 const RgbColor &color = RgbColor(1,1,1));
+
   virtual void 
   ExportBspSplitPlanes(const BspTree &tree);
@@ -101,13 +111,4 @@
   ExportLeavesGeometry(const BspTree &tree, const vector<BspLeaf *> &leaves);
 
-  bool
-  ExportRays(const RayContainer &rays,
-						 const float length=1000,
-						 const RgbColor &color = RgbColor(1,1,1));
-
-	bool
-	ExportRays(const VssRayContainer &rays,
-						 const RgbColor &color = RgbColor(1,1,1));
-		
   virtual void
   ExportBspViewCellPartition(const BspTree &tree, const int maxPvs = 0);
@@ -116,5 +117,13 @@
   ExportBspLeaves(const BspTree &tree, const int maxPvs = 0);
 
+
+  virtual void 
+  ExportBspSplits(const VspBspTree &tree, const bool exportDepth);
+
+  virtual void
+  ExportBspViewCellPartition(const VspBspTree &tree, const int maxPvs = 0);
+
 protected:
+
 
   virtual void
@@ -126,9 +135,13 @@
   bool
   ExportBspTreeRayDensity(const BspTree &tree);
+  
+  void
+  AddBoxToMesh(const AxisAlignedBox3 &box,
+			   Mesh *mesh);
 
-	void
-	AddBoxToMesh(const AxisAlignedBox3 &box,
-							 Mesh *mesh);
-		
+  void ExportBspNodeSplits(BspNode *root,
+						   const AxisAlignedBox3 &box, 
+						   const bool exportDepth,
+						   const bool epsilon);
 };
  
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h	(revision 446)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h	(revision 448)
@@ -135,4 +135,6 @@
   Vector3 GetNormalizedDir() const { return (mTermination - mOrigin)*mInvSize; }
 
+  float Length() const { return Distance(mOrigin, mTermination); }
+
   Vector3 Extrap(const float t) const {
 	return GetOrigin() + t * GetDir();
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp	(revision 447)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp	(revision 448)
@@ -9,5 +9,5 @@
 #include "VssRay.h"
 #include "RssTree.h"
-
+#include "ViewCellBsp.h"
 
 static bool useViewSpaceBox = true;//true;
@@ -474,21 +474,4 @@
   }
   
-  //-- construct BSP view cells
-  if (ViewCell::sHierarchy == ViewCell::BSP) 
-  {
-	const int bspSamples = min((int)mVssRays.size(), mBspConstructionSamples);
-	
-	Debug << "bpssamples: " << bspSamples << endl;
-	for (int i = 0; i < bspSamples; ++ i)
-	  bspRays.push_back(new Ray(*mVssRays[i]));
-	
-	//-- construct BSP tree using the samples
-	mBspTree = new BspTree(&mUnbounded);	
-	
-	ObjectContainer objects;
-	mSceneGraph->CollectObjects(&objects);
-	mBspTree->Construct(objects, bspRays);
-  }
-  
   rssTree = new RssTree;
   // viewcells = Construct(mVssRays);
@@ -567,13 +550,4 @@
 	}
 
-	// cast rays into BSP tree
-	if (ViewCell::sHierarchy == ViewCell::BSP) 
-	  {
-		for (int i = 0; i < (int)vssRays.size(); ++ i)
-		  {
-			CastRay(*mBspTree, *vssRays[i]);
-		  }
-	  }
-
 	samples+=num;
 	float pvs = rssTree->GetAvgPvsSize();
@@ -598,19 +572,4 @@
   delete rssTree;
   
-  if (ViewCell::sHierarchy == ViewCell::BSP) 
-	{
-	  Debug << mBspTree->GetStatistics();
-
-	  ObjectContainer objects;
-	  ExportSplits(objects, bspRays, 10000);
-	  ExportBspPvs(objects, bspRays, 10000);
-
-	  BspViewCellsStatistics stat;
-	  mBspTree->EvaluateViewCellsStats(stat);
-	  Debug << "original view cell partition:\n" << stat << endl;
-
-	  // clear BSP samples
-	  CLEAR_CONTAINER(bspRays);
-	}
 
   return true;
@@ -646,7 +605,7 @@
 	  BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
 	  // if ray not in unbounded space
-	  if (leaf->GetViewCell() != &mUnbounded)
+/*	  if (leaf->GetViewCell() != &mUnbounded)
 		contributingSamples += 
-		  leaf->GetViewCell()->GetPvs().AddSample(obj);
+		  leaf->GetViewCell()->GetPvs().AddSample(obj);*/
 	}
  
@@ -657,7 +616,7 @@
 		BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
 		
-		if (leaf->GetViewCell() != &mUnbounded)
+	/*	if (leaf->GetViewCell() != &mUnbounded)
 		  leaf->GetViewCell()->
-			AddPassingRay(ray, contributingSamples ? 1 : 0);
+			AddPassingRay(ray, contributingSamples ? 1 : 0);*/
 	  }
   
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp	(revision 446)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp	(revision 448)
@@ -384,6 +384,4 @@
   VssTree *vssTree = NULL;
 
-  RayContainer bspRays;
-
   while (totalSamples < mInitialSamples) {
 	int passContributingSamples = 0;
@@ -454,5 +452,6 @@
   mSceneGraph->CollectObjects(&objects);
 
-  //mViewCellsManager->Construct(mVssRays, objects);
+  // construct view cells
+  mViewCellsManager->Construct(objects, mVssRays, mViewSpaceBox);
 
   vssTree = new VssTree;
@@ -489,4 +488,8 @@
   int samples = 0;
   int pass = 0;
+
+  /// Rays used for post processing and visualizations.
+  RayContainer storedRays;
+  // cast view cell samples
   while (1) {
 	int num = mVssSamplesPerPass;
@@ -503,6 +506,5 @@
 	  num = GenerateImportanceRays(vssTree, num, rays);
 	}
-		
-		
+			
 	for (int i=0; i < rays.size(); i++)
 	  CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays);
@@ -528,13 +530,34 @@
 	}
 
-	// cast rays into BSP tree
-	/*if (ViewCell::sHierarchy == ViewCell::BSP) 
-	  {
-		for (int i = 0; i < (int)vssRays.size(); ++ i)
-		  {
-			CastRay(*mBspTree, *vssRays[i]);
-		  }
-	  }
-*/
+	//-- prepare traversal rays for view cell intersections
+	RayContainer passRays;
+
+	VssRayContainer::const_iterator it, it_end = vssRays.end();
+	
+	for (it = vssRays.begin(); it != it_end; ++ it)
+		passRays.push_back(new Ray(*(*it)));
+	
+	int sampleContributions = 0;
+	int contributingSamples = 0;
+
+	/// compute view cell contribution of rays
+	mViewCellsManager->ComputeSampleContributions(passRays,
+												  sampleContributions, 
+												  contributingSamples);
+	
+	//-- save rays for post processing
+	if (((int)storedRays.size() < mViewCellsManager->GetPostProcessSamples()) ||
+		((int)storedRays.size() < mViewCellsManager->GetVisualizationSamples()))
+	{
+		RayContainer::const_iterator it, it_end = passRays.end();
+
+		for (it = passRays.begin(); it != it_end; ++ it)
+			storedRays.push_back(new Ray(*(*it)));
+	}
+	else
+	{
+		CLEAR_CONTAINER(passRays);
+	}
+
 	samples+=num;
 	float pvs = vssTree->GetAvgPvsSize();
@@ -547,5 +570,11 @@
 	pass++;
   }
-	
+
+  //-- post process view cells
+  mViewCellsManager->PostProcess(objects, storedRays);
+  mViewCellsManager->Visualize(objects, storedRays);
+
+  CLEAR_CONTAINER(storedRays);
+
   delete vssTree;
 
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp	(revision 446)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp	(revision 448)
@@ -58,7 +58,4 @@
 	p->KdTreeStatistics(cout);
 	
-	// parse view cell hierarchy options
-	p->ParseViewCellsOptions();
-  
 	// parse view cells related options
 	p->PrepareViewCells();
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RayInfo.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RayInfo.cpp	(revision 437)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RayInfo.cpp	(revision 448)
@@ -45,15 +45,15 @@
 float RayInfo::ExtrapOrigin(const int axis) const 
 {
-	return mRay->GetOrigin(axis) + GetMinT()*mRay->GetDir(axis);
+	return mRay->GetOrigin(axis) + GetMinT() * mRay->GetDir(axis);
 }
 		
 float RayInfo::ExtrapTermination(const int axis) const 
 {
-	return mRay->GetOrigin(axis) + GetMaxT()*mRay->GetDir(axis);
+	return mRay->GetOrigin(axis) + GetMaxT() * mRay->GetDir(axis);
 }
 
 Vector3 RayInfo::ExtrapOrigin() const 
 {
-	return mRay->GetOrigin() * mRay->GetDir();
+	return mRay->GetOrigin() + GetMinT() * mRay->GetDir();
 }
 		
@@ -134,4 +134,6 @@
 	t = splitPlane.FindT(mRay->GetOrigin(), mRay->GetTermination());
 
+//	Debug << "t: " << t << " mint " << GetMinT() << " maxt " << GetMaxT() << " orig: " << mRay->GetOrigin() <<
+//		" term " << mRay->GetTermination() << endl;
 	if ((t < GetMinT()) || (t > GetMaxT()))
 		return splitPlane.Side(ExtrapOrigin());
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp	(revision 436)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp	(revision 448)
@@ -42,9 +42,10 @@
 void Polygon3::Split(const Plane3 &partition, 
 					 Polygon3 &front, 
-					 Polygon3 &back)
+					 Polygon3 &back,
+					 const float epsilon)
 {
 	Vector3 ptA = mVertices.back();
 	
-	int sideA = partition.Side(ptA, Vector3::sDistTolerance);
+	int sideA = partition.Side(ptA, epsilon);
 	
 	VertexContainer::const_iterator it;
@@ -53,9 +54,10 @@
 
 	bool foundSplit = false;
-	// find line - plane intersections
+
+	//-- find line - plane intersections
 	for (it = mVertices.begin(); it != mVertices.end(); ++ it)
 	{
 		Vector3 ptB = *it;
-		int sideB = partition.Side(ptB, Vector3::sDistTolerance);
+		int sideB = partition.Side(ptB, epsilon);
 	
 		// vertices on different sides => split
@@ -68,6 +70,5 @@
 			
 				// test if split point not too close to previous split point
-				if (!foundSplit || 
-					(SqrDistance(splitPt, lastSplit) > Vector3::sDistToleranceSqrt))
+				if (!foundSplit || (!EpsilonEqualV3(splitPt, lastSplit, epsilon)))
 				{
 					// add vertex to both polygons
@@ -89,6 +90,5 @@
 				// test if split point not too close to other split point
 					// test if split point not too close to previous split point
-				if (!foundSplit || 
-					(SqrDistance(splitPt, lastSplit) > Vector3::sDistToleranceSqrt))
+				if (!foundSplit || (!EpsilonEqualV3(splitPt, lastSplit, epsilon)))
 				{
 					// add vertex to both polygons
@@ -125,7 +125,8 @@
 }
 
-int Polygon3::Side(const Plane3 &plane) const
-{
-	int classification = ClassifyPlane(plane);
+int Polygon3::Side(const Plane3 &plane, 
+				   const float epsilon) const
+{
+	int classification = ClassifyPlane(plane, epsilon);
 	
 	if (classification == BACK_SIDE)
@@ -137,5 +138,6 @@
 }
 
-int Polygon3::ClassifyPlane(const Plane3 &plane, const float tolerance) const
+int Polygon3::ClassifyPlane(const Plane3 &plane, 
+							const float epsilon) const
 {
 	VertexContainer::const_iterator it;
@@ -148,5 +150,5 @@
 	for (it = mVertices.begin(); it != mVertices.end(); ++ it)
 	{
-		const int side = plane.Side(*it, tolerance);
+		const int side = plane.Side(*it, epsilon);
 		
 		if (side > 0)
@@ -179,8 +181,8 @@
 }
 
-int Polygon3::ClassifyPlane(const Plane3 &plane) const
-{
-	return ClassifyPlane(plane, Vector3::sDistTolerance);
-}
+/*int Polygon3::ClassifyPlane(const Plane3 &plane) const
+{
+	return ClassifyPlane(plane, Limits::Small);
+}*/
 
 
@@ -207,5 +209,5 @@
 }
 
-bool Polygon3::Valid() const
+bool Polygon3::Valid(const float epsilon) const
 {
 	if (mVertices.size() < 3)
@@ -225,7 +227,7 @@
 	for (it = mVertices.begin(); it != it_end; ++it)
 	{
-		if (!(SqrDistance(vtx, *it) > sDistToleranceSqrt))
-		{
-			Debug << "Malformed vertices:\n" << *this << endl;
+		if (!(EpsilonEqual(vtx, *it)))
+		{
+			//Debug << "Malformed vertices:\n" << *this << endl;
 			return false;
 		}
@@ -385,5 +387,7 @@
 }
 
-int Polygon3::ClassifyPlane(const PolygonContainer &polys, const Plane3 &plane)
+int Polygon3::ClassifyPlane(const PolygonContainer &polys, 
+							const Plane3 &plane,
+							const float epsilon)
 {
 	PolygonContainer::const_iterator it;
@@ -395,5 +399,5 @@
 	for (it = polys.begin(); it != polys.end(); ++ it)
 	{
-        const int cf = (*it)->ClassifyPlane(plane);
+        const int cf = (*it)->ClassifyPlane(plane, epsilon);
 		
         if (cf == FRONT_SIDE)
@@ -445,5 +449,10 @@
 
 	for (pit = cell.begin(); pit != cell.end(); ++ pit)
+	{
+		if ((*pit)->mVertices.size() < 3)
+			Debug << "ERROR" << endl;
+
 		area += (*pit)->GetArea();
+	}
 
 	return area;
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 444)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 448)
@@ -75,7 +75,7 @@
 
 	mViewCellsManager->ComputeSampleContributions(rays, 
-												  castRaysToViewCells,
 												  sampleContributions, 
-												  contributingSamples);
+												  contributingSamples,
+												  castRaysToViewCells);
 }
 
@@ -397,6 +397,6 @@
 	//-- several visualizations and statistics
 		 
-	//-- render simulation
-	cout << "\nevaluating bsp view cells render time before merge ... ";
+	//-- render simulation after merge
+	cout << "\nevaluating bsp view cells render time after merge ... ";
 		 
 	const SimulationStatistics ss = mViewCellsManager->SimulateRendering();
@@ -405,5 +405,6 @@
 	cout << ss << endl;
 	Debug << ss << endl;
-		 
+	
+
 	mViewCellsManager->Visualize(objects, mSampleRays);	
     
@@ -424,5 +425,6 @@
 		RayContainer::const_iterator it, it_end = newRays.end();
 
-		if ((int)mVssSampleRays.size() < mViewCellsManager->GetConstructionSamples())
+		if ((int)mVssSampleRays.size() < 
+			mViewCellsManager->GetConstructionSamples())
 		{
 			for (it = newRays.begin(); it != it_end; ++ it)
@@ -433,14 +435,12 @@
 		else
 		{
-			// construct view cells using the collected samples
-			cout << "building view cells from " << (int)mVssSampleRays.size() << " samples " << endl;
-
+			// construct view cells
 			mViewCellsManager->Construct(objects, mVssSampleRays);
-			
+		
 			// throw away samples 
-			CLEAR_CONTAINER(mVssSampleRays);
-		}
-	}
-	// Need (ordered) rays for post processing => collect new rays
+			//CLEAR_CONTAINER(mVssSampleRays);
+		}
+	}
+	// Need rays (with ordered intersections) for post processing => collect new rays
 	else if (((int)mSampleRays.size() < mViewCellsManager->GetPostProcessSamples()) ||
 			 ((int)mSampleRays.size() < mViewCellsManager->GetVisualizationSamples()))
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp	(revision 372)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp	(revision 448)
@@ -1662,5 +1662,5 @@
 
 // TODO: use a table to avoid normal and distance computations
-Polygon3 *AxisAlignedBox3::CrossSection(const Plane3 &plane)
+Polygon3 *AxisAlignedBox3::CrossSection(const Plane3 &plane) const
 {
 	Polygon3 *planePoly = new Polygon3();
@@ -1746,2 +1746,27 @@
 	return planePoly;
 }
+
+bool AxisAlignedBox3::GetRaySegment(const Ray &ray, 
+									float &minT, 
+									float &maxT) const
+{
+	maxT = 1e6;
+	minT = 0;
+	
+	// test with  bounding box
+	if (!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;
+}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.cpp	(revision 444)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.cpp	(revision 448)
@@ -68,7 +68,7 @@
 
 void ViewCellsManager::ComputeSampleContributions(const RayContainer &rays, 
-												  const bool castRays,
 												  int &sampleContributions,
-												  int &contributingSamples)
+												  int &contributingSamples,
+												  const bool castRays)
 {
 	// view cells not yet constructed
@@ -168,5 +168,5 @@
 SimulationStatistics ViewCellsManager::SimulateRendering() const
 {
-	if (!ViewCellsConstructed())
+	if (!ViewCellsConstructed() || !mRenderSimulator)
 		return SimulationStatistics();
 
@@ -240,6 +240,4 @@
 		return 0;
 
-	Debug << "Constructing bsp view cells" << endl;
-
 	int sampleContributions = 0;
 	
@@ -247,4 +245,8 @@
 
 	int limit = min (mConstructionSamples, (int)rays.size());
+
+	// construct view cells using the collected samples
+	Debug << "construcing bsp view cells from " 
+		  << limit << " samples " << endl;
 
     for (int i = 0; i < limit; ++ i)
@@ -332,9 +334,9 @@
 	int pvsSize = 0;
 
-	BspViewCellsStatistics stat;
-	mBspTree->EvaluateViewCellsStats(stat);
-	Debug << "original view cell partition:\n" << stat << endl;
-		
-	if (1) // export view cells
+	BspViewCellsStatistics vcStats;
+	mBspTree->EvaluateViewCellsStats(vcStats);
+	Debug << "original view cell partition:\n" << vcStats << endl;
+		
+	if (1) // export view cells before merge
 	{
 		cout << "exporting initial view cells (=leaves) ... ";
@@ -344,5 +346,5 @@
 		{
 			exporter->SetWireframe();
-			exporter->ExportBspLeaves(*mBspTree, stat.maxPvs);
+			exporter->ExportBspLeaves(*mBspTree, vcStats.maxPvs);
 			//exporter->ExportBspViewCellPartition(*mBspTree, 0);
 				
@@ -360,10 +362,19 @@
 		}
 		cout << "finished" << endl;
-	}
+
+		//-- render simulation
+		cout << "\nevaluating bsp view cells render time before merge ... ";
+		 
+		const SimulationStatistics ss = SimulateRendering();
+			 
+		cout << " finished" << endl;
+		cout << ss << endl;
+		Debug << ss << endl;
+	}
+
 
 	cout << "starting post processing using " << mPostProcessSamples << " samples ... ";
-		
+
 	long startTime = GetTime();
-
 
 	//-- merge or subdivide view cells
@@ -408,4 +419,9 @@
 	cout << "merged " << merged << " view cells in "
 		 << TimeDiff(startTime, GetTime()) *1e-3 << " secs" << endl;
+
+	// statistics after post process
+	vcStats.Reset();
+	mBspTree->EvaluateViewCellsStats(vcStats);
+	Debug << "post processed view cell partition:\n" << vcStats << endl;
 
 	return merged;
@@ -956,5 +972,5 @@
 mVspKdTree(vspKdTree)
 {
-   mRenderSimulator = NULL; // TODO: new VspKdRenderSimulator(vspKdTree);
+   mRenderSimulator = NULL;//new VspKdRenderSimulator(vspKdTree);
    InitRenderSimulator();
 }
@@ -983,23 +999,9 @@
 		return 0;
 
-	//if (castRay) // TODO
+	//if (castRay)
+		//mVspKdTree->CastRay(ray);
 
 	int sampleContributions = 0;
 
-/*
-    for (int i = 0; i < (int)rays.size(); ++ i)
-	{
-		Ray *ray = rays[i];
-
-		mBspTree->CastRay(*ray);
-			
-		Intersectable *term = 
-			!ray->intersections.empty() ? ray->intersections[0].mObject : NULL
-			;
-		sampleContributions += 
-			AddSampleContributions(*ray, ray->sourceObject.mObject, term);
-			
-	}
-*/
 	return sampleContributions;
 }
@@ -1121,6 +1123,6 @@
 mVspBspTree(vspBspTree)
 {
-   mRenderSimulator = NULL;// TODO new VspBspRenderSimulator(bspTree);
-   //InitRenderSimulator();
+	mRenderSimulator = new VspBspRenderSimulator(vspBspTree);
+	InitRenderSimulator();
 }
 
@@ -1182,5 +1184,5 @@
 		for (int j = 0; j < (int)ray.bspIntersections.size(); ++ j)
 		{ 	
-			VspBspLeaf *leaf = NULL;  //TODO ray.bspIntersections[j].mLeaf;
+			BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
 		
 			// if ray not in unbounded space
@@ -1206,5 +1208,5 @@
 		for (int j = 1; j < ((int)ray.bspIntersections.size() - 1); ++ j) 
 		{
-			VspBspLeaf *leaf = NULL;// TODO = ray.bspIntersections[j].mLeaf;
+			BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
 
 			if (leaf->GetViewCell() != mVspBspTree->GetRootCell())
@@ -1219,5 +1221,4 @@
 										const RayContainer &rays)
 {
-	return 0;// TODO
 	if (!ViewCellsConstructed())
 	{
@@ -1231,5 +1232,5 @@
 	int pvsSize = 0;
 
-	VspBspViewCellsStatistics vcStats;
+	BspViewCellsStatistics vcStats;
 	mVspBspTree->EvaluateViewCellsStats(vcStats);
 	Debug << "original view cell partition:\n" << vcStats << endl;
@@ -1243,5 +1244,5 @@
 		{
 			exporter->SetWireframe();
-			//exporter->ExportBspLeaves(*mVspBspTree, stat.maxPvs);
+			exporter->ExportBspViewCellPartition(*mVspBspTree, vcStats.maxPvs);
 				
 			if (0)
@@ -1284,10 +1285,10 @@
 		iit = ray->bspIntersections.begin();
 
-		VspBspLeaf *previousLeaf = NULL;//TODO(*iit).mLeaf;
+		BspLeaf *previousLeaf = (*iit).mLeaf;
 		++ iit;
 		
 		for (; iit != ray->bspIntersections.end(); ++ iit)
 		{
-			VspBspLeaf *leaf = NULL;//TODO(*iit).mLeaf;
+			BspLeaf *leaf = (*iit).mLeaf;
 
 			if (ShouldMerge(leaf, previousLeaf))
@@ -1307,4 +1308,9 @@
 		 << TimeDiff(startTime, GetTime()) *1e-3 << " secs" << endl;
 
+	// statistics after post process
+	vcStats.Reset();
+	mVspBspTree->EvaluateViewCellsStats(vcStats);
+	Debug << "post processed view cell partition:\n" << vcStats << endl;
+
 	return merged;
 }
@@ -1319,10 +1325,9 @@
 									   const RayContainer &sampleRays)
 {
-	return; // TODO
 	if (!ViewCellsConstructed())
 		return;
 
 	//-- recount pvs
-	VspBspViewCellsStatistics vcStats;
+	BspViewCellsStatistics vcStats;
 	mVspBspTree->EvaluateViewCellsStats(vcStats);
 
@@ -1330,5 +1335,5 @@
 	{
 		cout << "exporting view cells after merge ... ";
-		/*Exporter *exporter = Exporter::GetExporter("merged_view_cells.x3d");
+		Exporter *exporter = Exporter::GetExporter("merged_view_cells.x3d");
 
 		if (exporter)
@@ -1337,5 +1342,5 @@
 			delete exporter;
 		}
-		*/
+		
 		cout << "finished" << endl;
 	}	
@@ -1343,5 +1348,5 @@
 	//-- visualization of the BSP splits
 	bool exportSplits = false;
-		environment->GetBoolValue("BspTree.Visualization.exportSplits", exportSplits);
+	environment->GetBoolValue("BspTree.Visualization.exportSplits", exportSplits);
 	
 	if (exportSplits)
@@ -1352,5 +1357,5 @@
 	}
 
-//TODO	ExportVspBspPvs(objects, sampleRays);
+	ExportBspPvs(objects, sampleRays);
 }
 
@@ -1367,6 +1372,6 @@
 		exporter->SetForcedMaterial(m);
 		exporter->SetWireframe();
-		// TODO
-		//exporter->ExportBspSplits(*mVspBspTree, true);
+		
+		exporter->ExportBspSplits(*mVspBspTree, true);
 
 		// take forced material, else big scenes cannot be viewed
@@ -1412,10 +1417,9 @@
 	const int raysOut = min((int)sampleRays.size(), mVisualizationSamples);
 	cout << "visualization using " << mVisualizationSamples << " samples" << endl;
-	vector<Ray *> vcRays[leafOut];
-
-	if (0)
+	
+	if (1)
 	{
 		//-- some random view cells and rays for output
-		vector<VspBspLeaf *> vspBspLeaves;
+		vector<BspLeaf *> vspBspLeaves;
 
 		for (int i = 0; i < leafOut; ++ i)
@@ -1424,7 +1428,8 @@
 		for (int i = 0; i < (int)vspBspLeaves.size(); ++ i)
 		{
+			RayContainer vcRays;
 			cout << "creating output for view cell " << i << " ... ";
 			// check whether we can add the current ray to the output rays
-			for (int k = 0; k < raysOut; ++ k) 
+			/*for (int k = 0; k < raysOut; ++ k) 
 			{
 				Ray *ray = sampleRays[k];
@@ -1436,8 +1441,9 @@
 					if (vspBspLeaves[i]->GetViewCell() == leaf->GetViewCell()) 
 					{
-						vcRays[i].push_back(ray);
+						vcRays.push_back(ray);
 					}
 				}
 			}
+			*/
 
 			Intersectable::NewMail();
@@ -1474,14 +1480,15 @@
 
 			Debug << i << ": pvs size=" << (int)vc->GetPvs().GetSize() 
-					<< ", piercing rays=" << (int)vcRays[i].size() << endl;
+					<< ", piercing rays=" << (int)vcRays.size() << endl;
 
 			// export rays piercing this view cell
-			exporter->ExportRays(vcRays[i], 1000, RgbColor(0, 1, 0));
-
-			m.mDiffuseColor = RgbColor(1, 0, 0);
+			exporter->ExportRays(vcRays, 1000, RgbColor(1, 0, 0));
+			exporter->ExportRays(vspBspLeaves[i]->mVssRays, RgbColor(1, 1, 1));
+			
+			m.mDiffuseColor = RgbColor(0, 1, 1);
 			exporter->SetForcedMaterial(m);
 
-			// exporter->SetWireframe();
-			exporter->SetFilled();
+			exporter->SetWireframe();
+			//exporter->SetFilled();
 
 			// output PVS of view cell
@@ -1527,5 +1534,5 @@
 		{
 			cout << "creating output for view cell " << i << " ... ";
-			
+			RayContainer vcRays;
             Intersectable::NewMail();
 			BspViewCell *vc = dynamic_cast<BspViewCell *>(viewCells[i]);
@@ -1539,9 +1546,9 @@
 				for	(int j = 0; j < (int)ray->bspIntersections.size(); ++ j)
 				{
-					VspBspLeaf *leaf = NULL; // TODO = ray->bspIntersections[j].mLeaf;
+					BspLeaf *leaf = ray->bspIntersections[j].mLeaf;
 
 					if (vc == leaf->GetViewCell()) 
 					{
-						vcRays[i].push_back(ray);
+						vcRays.push_back(ray);
 					}
 				}
@@ -1572,9 +1579,9 @@
 			
 			Debug << i << ": pvs size=" << (int)vc->GetPvs().GetSize() 
-					<< ", piercing rays=" << (int)vcRays[i].size() << endl;
+					<< ", piercing rays=" << (int)vcRays.size() << endl;
 
 			
 			// export rays piercing this view cell
-			exporter->ExportRays(vcRays[i], 1000, RgbColor(0, 1, 0));
+			exporter->ExportRays(vcRays, 1000, RgbColor(0, 1, 0));
 	
 			m.mDiffuseColor = RgbColor(1, 0, 0);
@@ -1607,6 +1614,6 @@
 
 
-bool VspBspViewCellsManager::MergeVspBspLeafViewCells(VspBspLeaf *front, 
-													  VspBspLeaf *back) const
+bool VspBspViewCellsManager::MergeVspBspLeafViewCells(BspLeaf *front, 
+													  BspLeaf *back) const
 {
 	BspViewCell *viewCell = 
@@ -1623,9 +1630,8 @@
 	BspViewCell *bVc = back->GetViewCell();
 
-	// TODO
-	/*vector<VspBspLeaf *> fLeaves = fVc->mBspLeaves;
-	vector<VspBspLeaf *> bLeaves = bVc->mBspLeaves;
-
-	vector<VspBspLeaf *>::const_iterator it;
+	vector<BspLeaf *> fLeaves = fVc->mBspLeaves;
+	vector<BspLeaf *> bLeaves = bVc->mBspLeaves;
+
+	vector<BspLeaf *>::const_iterator it;
 	
 	for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
@@ -1642,9 +1648,9 @@
 	DEL_PTR(fVc);
 	DEL_PTR(bVc);
-*/
+
 	return true;
 }
 
-bool VspBspViewCellsManager::ShouldMerge(VspBspLeaf *front, VspBspLeaf *back) const
+bool VspBspViewCellsManager::ShouldMerge(BspLeaf *front, BspLeaf *back) const
 {
 	ViewCell *fvc = front->GetViewCell();
