Index: trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 344)
+++ trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env	(revision 349)
@@ -8,8 +8,8 @@
 #	filename glasgow1.x3d
 #	filename vienna.x3d
-	filename ../data/vienna/vienna-simple.x3d
+#	filename ../data/vienna/vienna-simple.x3d
 #	filename ../data/vienna/vienna-buildings.x3d
 #	filename ../data/vienna/viewcells-25-sel.x3d
-#	filename ../data/atlanta/atlanta2.x3d
+	filename ../data/atlanta/atlanta2.x3d
 #	filename ../data/soda/soda.dat
 #	filename ../data/soda/soda5.dat
@@ -56,6 +56,6 @@
 
 Sampling {
-	totalSamples	100
-	samplesPerPass  20
+	totalSamples	200000
+	samplesPerPass  2
 }
 
@@ -76,8 +76,8 @@
 BspTree {
 	Construction {
-	#	input fromRays
-		input fromViewCells
+		input fromRays
+	#	input fromViewCells
 	#	input fromSceneGeometry
-		samples 200000
+		samples 150000
 		sideTolerance 0.005
 	}
@@ -113,16 +113,17 @@
 	#splitPlaneStrategy 130
 	
-	splitPlaneStrategy 66
+	#splitPlaneStrategy 384
+	splitPlaneStrategy 130
 	
-	maxCandidates 50
+	maxCandidates 80
 	
 	Termination {
 		# autopartition
-		maxRays -1
+		maxRays 10
 		maxPolygons 0
 		maxDepth 100
 		
 		#axis aligned split
-		maxPolysForAxisAligned 100
+		maxPolysForAxisAligned 200
 		maxCostRatio 0.9
 		ct_div_ci 0.5
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h	(revision 333)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h	(revision 349)
@@ -23,10 +23,16 @@
 {
 public:
+	enum {BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT};
+
+	/** Default constructor creating an empty polygon.
+	*/
 	Polygon3();
-
+	/** Constructor creating a polygon from the vertices.
+	*/
 	Polygon3(const VertexContainer &vertices);
+	/** Creates a polygon and stores pointer to parent mesh
+		instance.
+	*/
 	Polygon3(MeshInstance *parent);
-
-	~Polygon3();
 
 	// creates an "infinite" polygon from this plane
@@ -89,12 +95,4 @@
 	int CastRay(const Ray &ray, float &t, const float nearestT);
 
-	/** Returns piercing rays container.
-	*/
-	RayContainer *GetPiercingRays();
-
-	/** Adds a ray to the ray container.
-	*/
-	void AddPiercingRay(Ray *ray);
-
 	/// vertices are connected in counterclockwise order.
 	VertexContainer mVertices;
@@ -107,8 +105,12 @@
 
 	/// Rays piercing this polygon
-	RayContainer *mPiercingRays;
+	RayContainer mPiercingRays;
 
-	static float sSideTolerance;
-	static float sSideToleranceSqrt;
+	/** The piercing rays of the polygon are inherited by the child fragments
+		@parm front_piece the front fragment inheriting the front rays
+		@param back_piece the back fragment inheriting the back rays
+	*/
+	void InheritRays(Polygon3 &front_piece, 
+					 Polygon3 &back_piece) const;
 };
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 344)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h	(revision 349)
@@ -169,15 +169,18 @@
 					  PolygonContainer &backPolys, 
 					  PolygonContainer &coincident,
-					  bool storePolys = false);
-
-	/** Splits the rays into front and back rays according to split plane
+					  const bool storePolys = false);
+
+	/** Subdivides the rays into front and back rays according to the split plane.
 		@param rays contains the rays to be split. The rays are 
-			   distributed to front and back rays.
+			   distributed into front and back rays.
 		@param frontRays returns rays on the front side of the plane
 		@param backRays returns rays on the back side of the plane
-	*/
-	void SplitRays(RayContainer &rays, 
-				   RayContainer &frontRays, 
-				   RayContainer &backRays);
+		@param box the bounding box of the scene
+		@returns the number of splits
+	*/
+	int SplitRays(RayContainer &rays, 
+			      RayContainer &frontRays, 
+				  RayContainer &backRays,
+				  const AxisAlignedBox3 &box);
 
 	/** Stores polygon in node or discards them according to storePolys.
@@ -192,16 +195,6 @@
 	}
 
-protected:	
-
-	/** The piercing rays of the polygon are inherited by the child fragments
-		@param poly the parent polygon
-		@parm front_piece the front fragment inheriting the front rays
-		@param back_piece the back fragment inheriting the back rays
-	*/
-	void InheritRays(const Polygon3 &poly, 
-					 Polygon3 &front_piece, 
-					 Polygon3 &back_piece);
-
-	enum {BACK_RAY, FRONT_RAY, SPLIT_RAY};
+protected:
+
 	/// Splitting plane corresponding to this node
 	Plane3 mPlane;
@@ -474,6 +467,6 @@
 							   PolygonContainer &coincident,
 							   RayContainer &rays,
-							   RayContainer &backRays,
-							   RayContainer &frontRays);
+							   RayContainer &frontRays,
+							    RayContainer &backRays);
 
 	/** Filters polygons down the tree.
@@ -571,4 +564,9 @@
 							 vector<SortableEntry> &splitCandidates) const;
 
+	/** Selects an axis aligned split plane.
+		Returns true if split is valied
+	*/
+	bool SelectAxisAlignedPlane(Plane3 &plane, const PolygonContainer &polys) const;
+
 	/// Pointer to the root of the tree
 	BspNode *mRoot;
@@ -642,9 +640,18 @@
 		contribution for a balanced tree.
 	*/
-	static float sLeastSplitsTable[4];
+	static float sLeastPolySplitsTable[4];
 	/** Evaluates split plane classification with respect to the plane's 
 		contribution for a minimum number splits in the tree.
 	*/
-	static float sBalancedTreeTable[4];
+	static float sBalancedPolysTable[4];
+	/** Evaluates split plane classification with respect to the plane's 
+		contribution for a minimum number of ray splits.
+	*/
+	static float sLeastRaySplitsTable[5];
+	/** Evaluates split plane classification with respect to the plane's 
+		contribution for balanced rays.
+	*/
+	static float sBalancedRaysTable[5];
+
 };
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h	(revision 318)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h	(revision 349)
@@ -36,4 +36,8 @@
   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
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 340)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 349)
@@ -24,4 +24,7 @@
 
   enum { NO_INTERSECTION=0, INTERSECTION_OUT_OF_LIMITS, INTERSECTION };
+
+  /// if ray is on back (front) side of plane, or goes from the front (back) to the back (front)
+  enum {FRONT, BACK, BACK_FRONT, FRONT_BACK, COINCIDENT};
 
   struct Intersection {
@@ -190,5 +193,5 @@
   /** Classifies ray with respect to the plane.
   */
-  int ClassifyPlane(const Plane3 &plane, float minT, float maxT);
+  int ClassifyPlane(const Plane3 &plane, const float minT, const float maxT) const;
 
 private:
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h	(revision 340)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h	(revision 349)
@@ -143,5 +143,5 @@
   AxisAlignedBox3 GetFaceBox(const int faceIndex);
 
-  void GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
+  int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
 
 	virtual ostream &Describe(ostream &s) {
@@ -160,5 +160,5 @@
   }
 
-  void GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
+  int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
 
 
@@ -215,5 +215,5 @@
   virtual int Type() const { return TRANSFORMED_MESH_INSTANCE; }
 
-  void GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
+  int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
 
 private:
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h	(revision 313)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h	(revision 349)
@@ -53,4 +53,9 @@
 	     ) = 0;
   
+  virtual bool
+  ExportRays(const RayContainer &rays,
+	     const float length=1000,
+	     const RgbColor &color = RgbColor(1,1,1)
+	     ) = 0;
 
   virtual void 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 342)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp	(revision 349)
@@ -43,9 +43,18 @@
 	contribution for a minimum number splits in the tree.
 */
-float BspTree::sLeastSplitsTable[] = {0, 0, 1, 0};
+float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
 /** Evaluates split plane classification with respect to the plane's 
 	contribution for a balanced tree.
 */
-float BspTree::sBalancedTreeTable[] = {-1, 1, 0, 0};
+float BspTree::sBalancedPolysTable[] = {1, -1, 0, 0};
+
+/** Evaluates split plane classification with respect to the plane's 
+	contribution for a minimum number of ray splits.
+*/
+float BspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
+/** Evaluates split plane classification with respect to the plane's 
+	contribution for balanced rays.
+*/
+float BspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
 
 /****************************************************************/
@@ -158,8 +167,11 @@
 }
 
-void BspInterior::SplitRays(RayContainer &rays, 
-							RayContainer &frontRays, 
-							RayContainer &backRays)
-{
+int BspInterior::SplitRays(RayContainer &rays, 
+						   RayContainer &frontRays, 
+						   RayContainer &backRays,
+						   const AxisAlignedBox3 &box)
+{
+	int splits = 0;
+
 	while (!rays.empty())
 	{
@@ -172,9 +184,10 @@
 	
 		// test with tree bounding box
-/*		if (!mBox.GetMinMaxT(*ray, &minT, &maxT))
+		if (!box.GetMinMaxT(*ray, &minT, &maxT))
 			continue;
+
 		if (minT < 0) // start ray from origin
 			minT = 0;
-*/
+
 		// bound ray
 		if ((ray->GetType() == Ray::LOCAL_RAY) && 
@@ -185,28 +198,34 @@
 		}
 		
-		int classification = ray->ClassifyPlane(mPlane, minT, maxT);
-		
+		const int classification = ray->ClassifyPlane(mPlane, minT, maxT);
+		
+		ray->SetId(classification);
+
 		switch (classification)
 		{
-		case Plane3::COINCIDENT:
+		case Ray::COINCIDENT:
 			break;
-		case Plane3::BACK_SIDE:
-			ray->SetId(BACK_RAY);
+		case Ray::BACK:
 			backRays.push_back(ray);
 			break;
-		case Plane3::FRONT_SIDE:
-			ray->SetId(FRONT_RAY);
+		case Ray::FRONT:
 			frontRays.push_back(ray);
 			break;
-		case Plane3::SPLIT:
-			ray->SetId(SPLIT_RAY);
+		case Ray::FRONT_BACK:
+			ray->SetId(Ray::FRONT_BACK);
+			backRays.push_back(ray);
+			++ splits;			
+			break;
+		case Ray::BACK_FRONT:
 			frontRays.push_back(ray);
-			backRays.push_back(ray);
-			//extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &maxt);
+			++ splits;
 			break;
 		default:
+			Debug << "Should not come here" << endl;
 			break;
 		}
 	}
+
+	return splits;
 }
 
@@ -215,5 +234,5 @@
 							   PolygonContainer &backPolys, 
 							   PolygonContainer &coincident,
-							   bool storePolys)
+							   const bool storePolys)
 {
 	Polygon3 *splitPoly = NULL;
@@ -232,5 +251,5 @@
 
 		// classify polygon
-		int classification = poly->ClassifyPlane(mPlane);
+		const int classification = poly->ClassifyPlane(mPlane);
 
 		Polygon3 *front_piece = NULL;
@@ -241,14 +260,14 @@
 		switch (classification)
 		{
-			case Plane3::COINCIDENT:
+			case Polygon3::COINCIDENT:
 				coincident.push_back(poly);
 				break;			
-			case Plane3::FRONT_SIDE:	
+			case Polygon3::FRONT_SIDE:	
 				frontPolys.push_back(poly);
 				break;
-			case Plane3::BACK_SIDE:
+			case Polygon3::BACK_SIDE:
 				backPolys.push_back(poly);
 				break;
-			case Plane3::SPLIT:
+			case Polygon3::SPLIT:
 				front_piece = new Polygon3(poly->mParent);
 				back_piece = new Polygon3(poly->mParent);
@@ -262,8 +281,8 @@
 				++ splits; // increase number of splits
 
-				//-- inherit rays
-				if (poly->mPiercingRays)
-					InheritRays(*poly, *front_piece, *back_piece);
-
+				//-- inherit rays from parent polygon
+				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())
@@ -280,6 +299,5 @@
 				Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
 #endif
-				ProcessPolygon(&poly, storePolys);
-				
+				ProcessPolygon(&poly, storePolys);			
 				break;
 			default:
@@ -290,26 +308,4 @@
 
 	return splits;
-}
-
-void BspInterior::InheritRays(const Polygon3 &poly, 
-						      Polygon3 &front_piece, 
-						      Polygon3 &back_piece)
-{
-	RayContainer *rays = poly.mPiercingRays;
-
-	RayContainer::const_iterator it, it_end = rays->end();
-
-	for (it = rays->begin(); it != it_end; ++ it)
-	{
-		switch((*it)->GetId())
-		{
-		case BACK_RAY:
-			back_piece.GetPiercingRays()->push_back(*it);
-			break;
-		case FRONT_RAY:
-			front_piece.GetPiercingRays()->push_back(*it);
-			break;
-		}
-	}
 }
 
@@ -664,5 +660,5 @@
 			if (it != facePolyMap.end()) 
 			{
-				(*it).second->AddPiercingRay(ray);
+				(*it).second->mPiercingRays.push_back(ray);
 			} 
 			else
@@ -670,5 +666,5 @@
 				Polygon3 *poly = new Polygon3(face, obj->GetMesh());
 				polys->push_back(poly);
-				poly->AddPiercingRay(ray);
+				poly->mPiercingRays.push_back(ray);
 
 				facePolyMap[face] = poly;
@@ -842,6 +838,6 @@
 									PolygonContainer &coincident,
 									RayContainer &rays,
-									RayContainer &backRays,
-									RayContainer &frontRays)
+									RayContainer &frontRays,
+									RayContainer &backRays)
 {
 	mStat.nodes += 2;
@@ -855,6 +851,6 @@
 #endif
 	
-	// partition rays
-	interior->SplitRays(rays, frontRays, backRays);
+	// subdivide rays into front and back rays
+	interior->SplitRays(rays, frontRays, backRays, mBox);
 	
 	// split polygons with split plane
@@ -987,12 +983,61 @@
 }
 
+bool BspTree::SelectAxisAlignedPlane(Plane3 &plane, 
+									 const PolygonContainer &polys) const
+{
+	AxisAlignedBox3 box;
+	box.Initialize();
+	
+	// create bounding box of region
+	Polygon3::IncludeInBox(polys, box);
+        
+	int objectsBack = 0, objectsFront = 0;
+	int axis = 0;
+	float costRatio = MAX_FLOAT;
+	Vector3 position;
+
+	//-- area subdivision
+	for (int i = 0; i < 3; ++ i) 
+	{
+		float p = 0;
+		float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
+		
+		if (r < costRatio)
+		{
+			costRatio = r;
+			axis = i;
+			position = p;
+		}
+	}
+	
+	if (costRatio >= sMaxCostRatio)
+		return false;
+
+	Vector3 norm(0,0,0); norm[axis] = 1.0f;
+	plane = Plane3(norm, position);
+
+	return true;
+}
+
 Plane3 BspTree::SelectPlane(BspLeaf *leaf, 
 							PolygonContainer &polys, 
 							const RayContainer &rays)
 {
-	if (polys.size() == 0)
-	{
-		Debug << "Warning: No split candidate available\n";
-		return Plane3();
+	if (polys.empty())
+	{
+		Debug << "Warning: No autopartition polygon candidate available\n";
+	
+		// return axis aligned split
+		AxisAlignedBox3 box;
+		box.Initialize();
+	
+		// create bounding box of region
+		Polygon3::IncludeInBox(polys, box);
+
+		const int axis = box.Size().DrivingAxis();
+		const Vector3 position = (box.Min()[axis] + box.Max()[axis])*0.5f;
+
+		Vector3 norm(0,0,0); norm[axis] = 1.0f;
+		return Plane3(norm, position);
 	}
 
@@ -1000,33 +1045,7 @@
 		((int)polys.size() > sTermMaxPolysForAxisAligned))
 	{
-		AxisAlignedBox3 box;
-		box.Initialize();
-	
-		// todo: area subdivision
-		Polygon3::IncludeInBox(polys, box);
-        
-		int objectsBack = 0, objectsFront = 0;
-		int axis = 0;
-		float costRatio = MAX_FLOAT;
-		Vector3 position;
-
-		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 < sMaxCostRatio)
-		{
-			Vector3 norm(0,0,0); norm[axis] = 1.0f;
-			return Plane3(norm, position);
-		}
+		Plane3 plane;
+		if (SelectAxisAlignedPlane(plane, polys))
+			return plane;
 	}
 
@@ -1074,5 +1093,5 @@
 int BspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
 {
-	int candidateIdx = Random(currentIdx --);
+	const int candidateIdx = Random(currentIdx --);
 
 	// swap candidates to avoid testing same plane 2 times
@@ -1101,9 +1120,9 @@
 	// needed for balanced view cells criterium
 	ViewCell::NewMail();
-	int backId = ViewCell::sMailId;
+	const int backId = ViewCell::sMailId;
 	ViewCell::NewMail();
-	int frontId = ViewCell::sMailId;
+	const int frontId = ViewCell::sMailId;
 	ViewCell::NewMail();
-	int frontAndBackId = ViewCell::sMailId;
+	const int frontAndBackId = ViewCell::sMailId;
 
 	PolygonContainer::const_iterator it, it_end = polys.end();
@@ -1111,15 +1130,15 @@
 	for (it = polys.begin(); it != it_end; ++ it)
 	{
-		int classification = (*it)->ClassifyPlane(candidatePlane);
+		const int classification = (*it)->ClassifyPlane(candidatePlane);
 
 		if (sSplitPlaneStrategy & BALANCED_POLYS)
-			sumBalancedPolys += sBalancedTreeTable[classification];
+			sumBalancedPolys += sBalancedPolysTable[classification];
 		
 		if (sSplitPlaneStrategy & LEAST_SPLITS)
-			sumSplits += sLeastSplitsTable[classification];
+			sumSplits += sLeastPolySplitsTable[classification];
 
 		if (sSplitPlaneStrategy & LARGEST_POLY_AREA)
 		{
-			if (classification == Plane3::COINCIDENT)
+			if (classification == Polygon3::COINCIDENT)
 				sumPolyArea += (*it)->GetArea();
 			//totalArea += area;
@@ -1128,7 +1147,7 @@
 		if (sSplitPlaneStrategy & BLOCKED_RAYS)
 		{
-			float blockedRays = (float)(*it)->GetPiercingRays()->size();
-		
-			if (classification == Plane3::COINCIDENT)
+			const float blockedRays = (float)(*it)->mPiercingRays.size();
+		
+			if (classification == Polygon3::COINCIDENT)
 				sumBlockedRays += blockedRays;
 			
@@ -1143,5 +1162,5 @@
 			// assure that we only count a view cell 
 			// once for the front and once for the back side of the plane
-			if (classification == Plane3::FRONT_SIDE)
+			if (classification == Polygon3::FRONT_SIDE)
 			{
 				if ((viewCell->mMailbox != frontId) && 
@@ -1158,5 +1177,5 @@
 				}
 			}
-			else if (classification == Plane3::BACK_SIDE)
+			else if (classification == Polygon3::BACK_SIDE)
 			{
 				if ((viewCell->mMailbox != backId) &&
@@ -1190,5 +1209,5 @@
 	if (sSplitPlaneStrategy & BLOCKED_RAYS)
 		if (totalBlockedRays != 0)
-			val += sBlockedRaysFactor * sumBlockedRays / totalBlockedRays;
+			val += sBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;
 
 	if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
@@ -1218,8 +1237,8 @@
 	
 		// test with tree bounding box
-		/*if (!mBox.GetMinMaxT(*ray, &minT, &maxT))
+		if (!mBox.GetMinMaxT(*ray, &minT, &maxT))
 			continue;
 		if (minT < 0) // start ray from origin
-			minT = 0;*/
+			minT = 0;
 
 		// bound ray
@@ -1231,14 +1250,15 @@
 		}
 
-		int classification = ray->ClassifyPlane(candidatePlane, minT, maxT);
+		const int classification = 
+			ray->ClassifyPlane(candidatePlane, minT, maxT);
 
 		if (sSplitPlaneStrategy & LEAST_RAY_SPLITS)
 		{
-			sumBalancedRays += sBalancedTreeTable[classification];
+			sumBalancedRays += sBalancedRaysTable[classification];
 		}
 		
 		if (sSplitPlaneStrategy & BALANCED_RAYS)
 		{
-			sumRaySplits += sLeastSplitsTable[classification];
+			sumRaySplits += sLeastRaySplitsTable[classification];
 		}
 	}
@@ -1332,6 +1352,6 @@
 	environment->GetBoolValue("BspTree.storeSplitPolys", sStoreSplitPolys);
 
-	environment->GetFloatValue("BspTree.Construction.sideTolerance", Polygon3::sSideTolerance);
-	Polygon3::sSideToleranceSqrt = Polygon3::sSideTolerance * Polygon3::sSideTolerance;
+	environment->GetFloatValue("BspTree.Construction.sideTolerance", Vector3::sDistTolerance);
+	Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
 
     Debug << "BSP max depth: " << sTermMaxDepth << endl;
@@ -1379,4 +1399,5 @@
 		{
 			BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
 			nodeStack.push(interior->GetBack());
 			nodeStack.push(interior->GetFront());
@@ -1417,10 +1438,10 @@
 	}
 	
-#ifdef _DEBUG
+//#ifdef _DEBUG
 	Debug << "BSP stats: "
 		  << "Depth: " << data.mDepth << " (max: " << sTermMaxDepth << "), "
 		  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << sTermMaxPolygons << "), "
 		  << "#rays: " << (int)data.mRays->size() << " (max: " << sTermMaxRays << ")" << endl;
-#endif
+//#endif
 }
 
@@ -1520,4 +1541,8 @@
 			entp = extp;
 			mint = maxt; // NOTE: need this?
+
+			 if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
+				break;
+
 			BspRayTraversalData &s = tStack.top();
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp	(revision 340)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp	(revision 349)
@@ -2,4 +2,6 @@
 #include "Vector3.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/Plane3.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Plane3.h	(revision 327)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Plane3.h	(revision 349)
@@ -30,6 +30,4 @@
   }
 
-  // relation of objects with the plane
-  enum {BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT};
 
   float Distance(const Vector3 &v) const {
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 328)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 349)
@@ -158,26 +158,30 @@
 }
 
-int Ray::ClassifyPlane(const Plane3 &plane, float minT, float maxT)
+int Ray::ClassifyPlane(const Plane3 &plane, 
+					   const float minT, 
+					   const float maxT) const
 {
-	Vector3 entp = Extrap(minT);
-	Vector3 extp = Extrap(maxT);
+	const Vector3 entp = Extrap(minT);
+	const Vector3 extp = Extrap(maxT);
   
-	int entSide = plane.Side(entp);
-	int extSide = plane.Side(extp);
+	const int entSide = plane.Side(entp);
+	const int extSide = plane.Side(extp);
 
 	if ((entSide == 0) && (extSide == 0)) 
- 	{
-		return Plane3::COINCIDENT;
+	{
+		return COINCIDENT;
 	}
-	else if ((entSide <= 0) && (entSide <= 0))
+	else if ((entSide <= 0) && (extSide <= 0))
 	{
-		return Plane3::BACK_SIDE;
+		return BACK;
 	}
-	else if ((entSide >= 0) && (entSide >= 0))
+	else if ((entSide >= 0) && (extSide >= 0))
 	{
-		return Plane3::FRONT_SIDE;
+		return FRONT;
 	}
+	else if ((entSide <= 0) && (extSide >= 0))
+		return BACK_FRONT;
 	
-	return Plane3::SPLIT;
+	return FRONT_BACK;
 }
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h	(revision 329)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h	(revision 349)
@@ -76,4 +76,9 @@
   ExportBspSplits(const BspTree &tree);
 
+  bool
+  ExportRays(const RayContainer &rays,
+			 const float length=1000,
+			 const RgbColor &color = RgbColor(1,1,1));
+
   static ViewCellContainer foundViewCells; // todo: remove this
 protected:
Index: trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h	(revision 335)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h	(revision 349)
@@ -64,4 +64,6 @@
 
 	bool BuildBspTree();
+
+	void ExportSplits(const ObjectContainer &objects);
 };
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp	(revision 333)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp	(revision 349)
@@ -72,4 +72,54 @@
   stream<<"</IndexedLineSet>"<<endl;
   stream<<"</Shape>"<<endl;
+  return true;
+}
+
+bool
+X3dExporter::ExportRays(const RayContainer &rays,
+						const float length,
+						const RgbColor &color)
+{
+  RayContainer::const_iterator ri = rays.begin();
+
+  stream<<"<Shape>"<<endl;
+  stream<<"<Appearance>"<<endl;
+  stream<<"<Material ambientColor=\""<<color.r<<" "<<color.g<<" "<<color.b<<
+    "\" />"<<endl;
+  stream<<"</Appearance>"<<endl;
+  
+  stream<<"<IndexedLineSet coordIndex=\""<<endl;
+
+  int index = 0;
+  for (; ri != rays.end(); ri++) {
+    stream<<index<<" "<<index+1<<" -1\n";
+    index+=2;
+  }
+  
+  stream<<"\" >"<<endl;
+  
+  stream<<"<Coordinate  point=\""<<endl;
+  
+  ri = rays.begin();
+  for (; ri != rays.end(); ri++) {
+    Vector3 a = (*ri)->GetLoc();
+    
+    Vector3 b;
+    if (length < 0)
+      b = (*ri)->GetLoc() - length*(*ri)->GetDir();
+    else
+      if ((*ri)->intersections.size()==0)
+	b = (*ri)->GetLoc() + length*(*ri)->GetDir();
+      else
+	b = (*ri)->Extrap((*ri)->intersections[0].mT);
+    
+    stream<<a.x<<" "<<a.y<<" "<<a.z<<" ,";
+    stream<<b.x<<" "<<b.y<<" "<<b.z<<" ,\n";
+  }
+  
+  stream<<"\" >"<<endl;
+  stream<<"</Coordinate>"<<endl;
+  stream<<"</IndexedLineSet>"<<endl;
+  stream<<"</Shape>"<<endl;
+
   return true;
 }
@@ -668,5 +718,5 @@
 				Polygon3 *backPoly = new Polygon3();
 
-				if (planePoly->ClassifyPlane(*tData.mPlanes[i]) == Plane3::SPLIT)
+				if (planePoly->ClassifyPlane(*tData.mPlanes[i]) == Polygon3::SPLIT)
 				{
 					planePoly->Split(*tData.mPlanes[i], *frontPoly, *backPoly, splitPts);
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h	(revision 340)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h	(revision 349)
@@ -41,5 +41,5 @@
   virtual int Type() const = 0;
 
-  virtual void GetRandomSurfacePoint(Vector3 &point, Vector3 &normal) = 0;
+  virtual int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal) = 0;
 
 	virtual ostream &Describe(ostream &s) = 0;
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp	(revision 344)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp	(revision 349)
@@ -6,21 +6,18 @@
 #include "Ray.h"
 
-float Polygon3::sSideTolerance = 0.002f;
-float Polygon3::sSideToleranceSqrt = 0.000004f;
-
 Polygon3::Polygon3(): 
-mMaterial(NULL), mParent(NULL), mPiercingRays(NULL)
+mMaterial(NULL), mParent(NULL), mPiercingRays(0)
 {}
 
 Polygon3::Polygon3(const VertexContainer &vertices): 
-mVertices(vertices), mMaterial(NULL), mParent(NULL), mPiercingRays(NULL)
+mVertices(vertices), mMaterial(NULL), mParent(NULL), mPiercingRays(0)
 {}
 
 Polygon3::Polygon3(MeshInstance *parent): 
-mMaterial(NULL), mParent(parent), mPiercingRays(NULL)
+mMaterial(NULL), mParent(parent), mPiercingRays(0)
 {}
 
 Polygon3::Polygon3(Face *face, Mesh *parentMesh):
-mMaterial(NULL), mParent(NULL), mPiercingRays(NULL)
+mMaterial(NULL), mParent(NULL), mPiercingRays(0)
 {	
 	VertexIndexContainer::iterator it = face->mVertexIndices.begin();
@@ -30,9 +27,4 @@
 		mMaterial = parentMesh->mMaterial;
 	}
-}
-
-Polygon3::~Polygon3()
-{
-	DEL_PTR(mPiercingRays);
 }
 
@@ -55,5 +47,5 @@
 	Vector3 ptA = mVertices.back();
 	
-	int sideA = partition.Side(ptA, sSideTolerance);
+	int sideA = partition.Side(ptA, Vector3::sDistTolerance);
 	
 	VertexContainer::const_iterator it;
@@ -64,5 +56,5 @@
 	{
 		Vector3 ptB = *it;
-		int sideB = partition.Side(ptB, sSideTolerance);
+		int sideB = partition.Side(ptB, Vector3::sDistTolerance);
 	
 		// vertices on different sides => split
@@ -76,5 +68,5 @@
 				// test if split point not too close to previous split point
 				if (!foundSplit || 
-					(SqrDistance(splitPt, splitPts.back()) > sSideToleranceSqrt))
+					(SqrDistance(splitPt, splitPts.back()) > Vector3::sDistToleranceSqrt))
 				{
 					// add vertex to both polygons
@@ -97,5 +89,5 @@
 					// test if split point not too close to previous split point
 				if (!foundSplit || 
-					(SqrDistance(splitPt, splitPts.back()) > sSideToleranceSqrt))
+					(SqrDistance(splitPt, splitPts.back()) > Vector3::sDistToleranceSqrt))
 				{
 					// add vertex to both polygons
@@ -136,7 +128,7 @@
 	int classification = ClassifyPlane(plane);
 	
-	if (classification == Plane3::BACK_SIDE)
+	if (classification == BACK_SIDE)
 		return -1;
-	else if (classification == Plane3::FRONT_SIDE)
+	else if (classification == FRONT_SIDE)
 		return 1;
 
@@ -155,5 +147,5 @@
 	for (it = mVertices.begin(); it != mVertices.end(); ++ it)
 	{
-		int side = plane.Side(*it, sSideTolerance);
+		int side = plane.Side(*it, Vector3::sDistTolerance);
 		
 		if (side > 0)
@@ -165,10 +157,10 @@
 		if (onFrontSide && onBackSide) // split 
 		{
-			return Plane3::SPLIT;
+			return SPLIT;
 		}
 		// 3 vertices enough to decide coincident
 		else if (((++ count) >= 3) && !onFrontSide && !onBackSide)
-		{   //if (DotProd(plane.mNormal, GetNormal()) > 0)
-			return Plane3::COINCIDENT; 
+		{   
+			return COINCIDENT; 
 		}
 	}
@@ -176,12 +168,12 @@
 	if (onBackSide)
 	{
-		return Plane3::BACK_SIDE;
+		return BACK_SIDE;
 	}
 	else if (onFrontSide)
 	{
-		return Plane3::FRONT_SIDE;
-	}
-	//if (DotProd(plane.mNormal, GetNormal()) < 0) return Plane3::FRONT_SIDE;
-	return Plane3::COINCIDENT; // plane and polygon are coincident
+		return FRONT_SIDE;
+	}
+	
+	return COINCIDENT; // plane and polygon are coincident
 }
 
@@ -227,5 +219,5 @@
 	for (it = mVertices.begin(); it != it_end; ++it)
 	{
-		if (!(SqrDistance(vtx, *it) > SIDE_TOLERANCE_SQRD))
+		if (!(SqrDistance(vtx, *it) > sDistToleranceSqrt))
 		{
 			Debug << "Malformed vertices:\n" << *this << endl;
@@ -356,13 +348,29 @@
 }
 
-RayContainer *Polygon3::GetPiercingRays()
-{
-	if (!mPiercingRays)
-		mPiercingRays = new RayContainer();
-	return mPiercingRays;
-}
-
-void Polygon3::AddPiercingRay(Ray *ray)
-{
-	GetPiercingRays()->push_back(ray);
-}
+void Polygon3::InheritRays(Polygon3 &front_piece, 
+					       Polygon3 &back_piece) const
+{
+	RayContainer::const_iterator it, 
+		it_end = mPiercingRays.end();
+
+	for (it = mPiercingRays.begin(); it != it_end; ++ it)
+	{
+		switch((*it)->GetId())
+		{
+		case Ray::BACK:
+			back_piece.mPiercingRays.push_back(*it);
+			break;
+		case Ray::FRONT:
+			front_piece.mPiercingRays.push_back(*it);
+			break;
+		case Ray::FRONT_BACK:
+			back_piece.mPiercingRays.push_back(*it);
+			break;
+		case Ray::BACK_FRONT:
+			front_piece.mPiercingRays.push_back(*it);
+			break;
+		default:
+			break;
+		}
+	}
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 344)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 349)
@@ -55,7 +55,6 @@
 	DEL_PTR(mBspTree);
 	mBspTree = new BspTree(&mUnbounded);
-
 	ObjectContainer objects;
-	int prunedCells = 0;
+	
 	switch (BspTree::sConstructionMethod)
 	{
@@ -76,9 +75,6 @@
 		mBspTree->SetGenerateViewCells(true);
 		mBspTree->Construct(mSampleRays);
-//		prunedCells = mBspTree->MergeViewCells();
 		
 		mBspTree->CollectViewCells(mViewCells);
-
-//		Debug << "Collected " << (int)mViewCells.size() << " cells (" << prunedCells << " pruned)" << endl;
 		break;
 	default:
@@ -169,19 +165,19 @@
 
 	if (mViewCellsType == BSP_VIEW_CELLS)
+	{
+		// cast ray to BSP tree to get intersection with view cells
+		if (mBspTree)
 		{
-			// cast ray to BSP tree to get intersection with view cells
-			if (mBspTree)
-				{
-					mBspTree->CastRay(ray);
+			mBspTree->CastRay(ray);
 					
-					sampleContributions += AddObjectSamples(object, ray);
+			sampleContributions += AddObjectSamples(object, ray);
 				
-					if (!ray.intersections.empty()) // second intersection found
-						{
-							sampleContributions += 
-								AddObjectSamples(ray.intersections[0].mObject, ray);
-						}
-				}
+			if (!ray.intersections.empty()) // second intersection found
+			{
+				sampleContributions += 
+					AddObjectSamples(ray.intersections[0].mObject, ray);
+			}
 		}
+	}
 	else
 		{
@@ -319,4 +315,9 @@
 		int maxTimeIdx = 0;
 
+		bool collectSamplesForBsp = 
+			(mViewCellsType == BSP_VIEW_CELLS) &&
+			(BspTree::sConstructionMethod == BspTree::FROM_RAYS) &&
+			(totalSamples < mBspConstructionSamples);
+			
 		cout << "totalSamples: "  << totalSamples << endl;
 
@@ -382,5 +383,5 @@
 			}
 			
-			object->GetRandomSurfacePoint(point, normal);
+			int faceIndex = object->GetRandomSurfacePoint(point, normal);
 			
 			long samplesPerObjStart = GetTime();
@@ -418,6 +419,14 @@
 					
 					// construct a ray
-					SetupRay(ray, point, direction, Ray::LOCAL_RAY);
-					
+					SetupRay(ray, point, direction, 
+						collectSamplesForBsp ? Ray::LINE_SEGMENT : Ray::LOCAL_RAY);
+			
+					// also add origin to sample in order to extract it as input polygons
+					if (0 && collectSamplesForBsp)
+					{
+						MeshInstance *meshInst = dynamic_cast<MeshInstance *>(object);
+						ray.intersections.push_back(Ray::Intersection(0.0, meshInst, faceIndex));
+					}
+
 					sampleContributions = CastRay(object, ray);
 
@@ -443,13 +452,12 @@
 					if (mViewCellsType == BSP_VIEW_CELLS) 
 					{
-                        // save rays for bsp tree construction
-						if ((BspTree::sConstructionMethod = BspTree::FROM_RAYS) &&
-							(totalSamples < mBspConstructionSamples))
+						// save rays for bsp tree construction
+						if (collectSamplesForBsp)
 						{
-							mSampleRays.push_back(new Ray(ray));
+                            mSampleRays.push_back(new Ray(ray));
 						}
 						else
 						{
-                            // construct BSP tree using the samples
+							// construct BSP tree using the samples
 							if (!mBspTree)
 							{
@@ -463,5 +471,5 @@
 								Export("vc_bsptree.x3d", false, false, true);
 							}
-													
+								
 							// some random view cells for output
 							if (pvsViewCells.empty())
@@ -477,5 +485,5 @@
 							}
 							else
-							{
+							{	
 								// check whether we can add the current ray to the rays
 								for (int k = 0; k < (int)ray.viewCells.size(); ++ k)
@@ -484,5 +492,5 @@
 											vcRays[j].push_back(ray);							
 							}
-						}				
+						}
 					}
 		
@@ -549,5 +557,5 @@
 			"#AvgRayContrib\n" << avgRayContrib << endl;
 	}
-
+	
 	if (mViewCellsType == KD_VIEW_CELLS)	
 		cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
@@ -559,5 +567,5 @@
     exporter->ExportKdTree(*mKdTree);
 
-	if (mViewCellsType == BSP_VIEW_CELLS)	
+	if (mBspTree && (mViewCellsType == BSP_VIEW_CELLS))
 		exporter->ExportBspTree(*mBspTree);
 
@@ -579,7 +587,8 @@
     delete exporter;
   }
+
   //-- several visualizations and statistics
   if (1) {
-  if (mViewCellsType == BSP_VIEW_CELLS)
+  if (mBspTree && (mViewCellsType == BSP_VIEW_CELLS))
   {
 	  bool exportSplits = false;
@@ -588,37 +597,6 @@
 	  // export the bsp splits
 	  if (exportSplits)
-	  {
-		  Exporter *exporter = Exporter::GetExporter("bsp_splits.x3d");
-
-		  if (exporter) 
-		  {	
-			  Material m; 
-			  m.mDiffuseColor = RgbColor(1, 0, 0);
-			  exporter->SetForcedMaterial(m);
-			  exporter->SetWireframe();
-			  exporter->ExportBspSplits(*mBspTree);
-
-			  // take forced material, else big scenes cannot be viewed
-			  m.mDiffuseColor = RgbColor(0, 1, 0);
-			  exporter->SetForcedMaterial(m);
-			  exporter->SetFilled();
-
-			  exporter->ResetForcedMaterial();
-
-			  // export scene geometry
-			  if (0)
-			  {
-				  Material m;//= RandomMaterial();
-				  m.mDiffuseColor = RgbColor(0, 0, 1);
-				  exporter->SetForcedMaterial(m);
-                 
-				  for (int j = 0; j < objects.size(); ++ j)
-					   exporter->ExportIntersectable(objects[j]);
-			
-				  delete exporter;
-			  }
-		  }
-	  }
-
+	  	  ExportSplits(objects);
+		 
 	  for (int j = 0; j < pvsViewCells.size(); ++ j)
 	  {
@@ -666,5 +644,5 @@
 		  }
 		
-		  // output rest of the sobjects
+		  // output rest of the objects
 		  if (0)
 		  {
@@ -727,2 +705,60 @@
   return true;
 }
+
+
+void SamplingPreprocessor::ExportSplits(const ObjectContainer &objects)
+{
+	Exporter *exporter = Exporter::GetExporter("bsp_splits.x3d");
+
+	if (exporter) 
+	{	
+		cout << "exporting splits ... ";
+
+		Material m; 
+		m.mDiffuseColor = RgbColor(1, 0, 0);
+		exporter->SetForcedMaterial(m);
+		exporter->SetWireframe();
+		exporter->ExportBspSplits(*mBspTree);
+
+		// take forced material, else big scenes cannot be viewed
+		m.mDiffuseColor = RgbColor(0, 1, 0);
+		exporter->SetForcedMaterial(m);
+		exporter->SetFilled();
+
+		exporter->ResetForcedMaterial();
+		
+		// export rays
+		if (0)
+		{
+			RayContainer outRays;
+		
+			for (int i = 0; i < mSampleRays.size(); ++ i)
+			{
+				// only rays piercing geometry
+				if (!mSampleRays[i]->intersections.empty())
+					outRays.push_back(mSampleRays[i]);
+			}
+			if (BspTree::sConstructionMethod == BspTree::FROM_RAYS)
+			{
+				// export rays 
+				exporter->ExportRays(outRays, 1000, RgbColor(1, 1, 0));
+			}
+		}
+
+		// export scene geometry
+		if (1)
+		{
+			Material m;//= RandomMaterial();
+			m.mDiffuseColor = RgbColor(0, 1, 0);
+			exporter->SetForcedMaterial(m);
+            exporter->SetWireframe();
+
+			for (int j = 0; j < objects.size(); ++ j)
+				exporter->ExportIntersectable(objects[j]);
+		}
+
+		delete exporter;
+
+		cout << "finished" << endl;
+	}
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp	(revision 340)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp	(revision 349)
@@ -266,6 +266,5 @@
 }
 
-
-void
+int
 Mesh::GetRandomSurfacePoint(Vector3 &point, Vector3 &normal)
 {
@@ -285,6 +284,7 @@
 	
 	normal = GetFacePlane(faceIndex).mNormal;
-}
-
+
+	return faceIndex;
+}
 
 int
@@ -308,16 +308,17 @@
 
 
-void
+int
 MeshInstance::GetRandomSurfacePoint(Vector3 &point, Vector3 &normal)
 {
-  mMesh->GetRandomSurfacePoint(point, normal);
-}
-
-void
+  return mMesh->GetRandomSurfacePoint(point, normal);
+}
+
+int
 TransformedMeshInstance::GetRandomSurfacePoint(Vector3 &point, Vector3 &normal)
 {
-  mMesh->GetRandomSurfacePoint(point, normal);
+  int index = mMesh->GetRandomSurfacePoint(point, normal);
   point = mWorldTransform*point;
   normal = TransformNormal(mWorldTransform, normal);
+  return index;
 }
 
