Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.h	(revision 447)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.h	(revision 447)
@@ -0,0 +1,855 @@
+// ================================================================
+// $Id$
+// ****************************************************************
+// 
+// Initial coding by
+/**
+   @author Jiri Bittner
+*/
+
+#ifndef __RSSTREE_H__
+#define __RSSTREE_H__
+
+// Standard headers
+#include <iomanip>
+#include <vector>
+#include <functional>
+#include <stack>
+
+
+// User headers
+#include "VssRay.h"
+#include "AxisAlignedBox3.h"
+
+
+#define USE_KDNODE_VECTORS 1
+#define _RSS_STATISTICS
+#define _RSS_TRAVERSAL_STATISTICS
+
+
+#include "Statistics.h"
+#include "Ray.h"
+
+#include "Containers.h"
+// --------------------------------------------------------------
+// Static statistics for kd-tree search
+// --------------------------------------------------------------
+class RssStatistics :
+  public StatisticsBase
+{
+public:  
+  // total number of nodes
+  int nodes;
+  // number of splits along each of the axes
+  int splits[7];
+  // totals number of rays
+  int rays;
+  // initial size of the pvs
+  int initialPvsSize;
+  // total number of query domains
+  int queryDomains;
+  // total number of ray references
+  int rayRefs;
+
+  // max depth nodes
+  int maxDepthNodes;
+  // max depth nodes
+  int minPvsNodes;
+  int minRaysNodes;
+  // max ray contribution nodes
+  int maxRayContribNodes;
+  // max depth nodes
+  int minSizeNodes;
+  int maxCostRatioNodes;
+
+  // max number of rays per node
+  int maxRayRefs;
+  // number of dynamically added ray refs
+  int addedRayRefs;
+  // number of dynamically removed ray refs
+  int removedRayRefs;
+  
+  // Constructor
+  RssStatistics() {
+    Reset();
+  }
+
+  int Nodes() const {return nodes;}
+  int Interior() const { return nodes/2; }
+  int Leaves() const { return (nodes/2) + 1; }
+
+  void Reset() {
+    nodes = 0;
+    for (int i=0; i<7; i++)
+      splits[i] = 0;
+    rays = queryDomains = 0;
+    rayRefs = 0;
+    maxDepthNodes = 0;
+    minPvsNodes = 0;
+    minRaysNodes = 0;
+    maxRayRefs = 0;
+    addedRayRefs = removedRayRefs = 0;
+	initialPvsSize = 0;
+	maxRayContribNodes = 0;
+	minSizeNodes = 0;
+	maxCostRatioNodes = 0;
+  }
+
+  
+  void
+  Print(ostream &app) const;
+
+  friend ostream &operator<<(ostream &s, const RssStatistics &stat) {
+    stat.Print(s);
+    return s;
+  }
+  
+};
+
+
+
+
+class RssTreeInterior;
+
+
+// --------------------------------------------------------------
+// KD-tree node - abstract class
+// --------------------------------------------------------------
+class RssTreeNode
+{
+public:
+
+  struct RayInfo
+  {
+	// pointer to the actual ray
+	VssRay *mRay;
+	
+	RayInfo():mRay(NULL) {}
+	
+	RayInfo(VssRay *r):mRay(r)
+	{}
+	
+	enum {
+	  SOURCE_RAY = 0,
+	  TERMINATION_RAY,
+	  PASSING_RAY,
+	  CONTAINED_RAY
+	};
+	
+	int GetRayClass() const {
+	  return SOURCE_RAY;
+	}
+	
+	friend bool operator<(const RayInfo &a, const RayInfo &b) {
+	  return a.mRay < b.mRay;
+	}
+	
+	Vector3 GetOrigin() const {
+	  return mRay->GetOrigin();
+	}
+
+	float ExtrapOrigin(const int axis) const {
+	  return mRay->GetOrigin(axis);
+	}
+	
+	float ExtrapTermination(const int axis) const {
+	  return mRay->GetTermination(axis);
+	}
+
+	Vector3 Extrap(const float t) const {
+	  return mRay->Extrap(t);
+	}
+	
+	int
+	ComputeRaySide(const int axis,
+				   const float position
+				   ) const {
+	  
+	  // only compare the side of the origin
+	  float rpos = mRay->GetOrigin(axis);
+	  return (rpos <= position) ? -1 : 1;
+	}
+	
+  };
+  
+  typedef vector<RayInfo> RayInfoContainer;
+	
+  enum { EInterior, ELeaf };
+
+  /////////////////////////////////
+  // The actual data goes here
+  
+  // link to the parent
+  RssTreeInterior *parent;
+
+  enum {SPLIT_X=0, SPLIT_Y, SPLIT_Z, SPLIT_DIRX, SPLIT_DIRY, SPLIT_DIRZ};
+  
+  // splitting axis
+  char axis;
+	
+  // depth
+  unsigned char depth;
+  
+  //  short depth;
+  //
+  /////////////////////////////////
+  
+  inline RssTreeNode(RssTreeInterior *p);
+
+  
+  virtual ~RssTreeNode() {};
+  virtual int Type() const  = 0;
+  
+
+  bool IsLeaf() const { return axis == -1; }
+  
+  virtual void Print(ostream &s) const = 0;
+
+  virtual int GetAccessTime() {
+    return 0x7FFFFFF;
+  }
+
+  
+	
+};
+
+// --------------------------------------------------------------
+// KD-tree node - interior node
+// --------------------------------------------------------------
+class RssTreeInterior :
+  public RssTreeNode
+{
+public:
+  // plane in local modelling coordinates
+  float position;
+
+  // pointers to children
+  RssTreeNode *back, *front;
+
+  // the bbox of the node
+  AxisAlignedBox3 bbox;
+
+  // the bbox of the node
+  AxisAlignedBox3 dirBBox;
+  
+  // data for caching
+  long accesses;
+  long lastAccessTime;
+  
+  RssTreeInterior(RssTreeInterior *p):RssTreeNode(p),
+									  back(NULL),
+									  front(NULL),
+									  accesses(0),
+									  lastAccessTime(-1)
+  { }
+
+  virtual int GetAccessTime() {
+    return lastAccessTime;
+  }
+
+  void SetupChildLinks(RssTreeNode *b, RssTreeNode *f) {
+    back = b;
+    front = f;
+    b->parent = f->parent = this;
+  }
+
+  void ReplaceChildLink(RssTreeNode *oldChild, RssTreeNode *newChild) {
+    if (back == oldChild)
+      back = newChild;
+    else
+      front = newChild;
+  }
+
+  virtual int Type() const  { return EInterior; }
+  
+  virtual ~RssTreeInterior() {
+    if (back)
+      delete back;
+    if (front)
+      delete front;
+  }
+  
+  virtual void Print(ostream &s) const {
+    if (axis == 0)
+      s<<"x ";
+    else
+      if (axis == 1)
+		s<<"y ";
+      else
+		s<<"z ";
+    s<<position<<" ";
+    back->Print(s);
+    front->Print(s);
+  }
+
+  
+	
+  int ComputeRaySide(const
+					 RayInfo &rayData
+					 ) {
+	return rayData.ComputeRaySide(axis, position);
+  }
+
+};
+
+
+// --------------------------------------------------------------
+// KD-tree node - leaf node
+// --------------------------------------------------------------
+class RssTreeLeaf :
+  public RssTreeNode
+{
+private:
+  int mPvsSize;
+public:
+  static int mailID;
+  int mailbox;
+  
+  RayInfoContainer rays;
+  int mPassingRays;
+	
+  bool mValidPvs;
+  float mEntropyImportance;
+  float mPvsEntropy;
+  float mRayLengthEntropy;
+  
+	
+  RssTreeLeaf(RssTreeInterior *p,
+			  const int nRays
+			  ):RssTreeNode(p), rays(), mPvsSize(0), mPassingRays(0), mValidPvs(false) {
+    rays.reserve(nRays);
+  }
+  
+  virtual ~RssTreeLeaf() { }
+
+  virtual int Type() const  { return ELeaf; }
+
+  virtual void Print(ostream &s) const {
+    s<<
+	  "L: rays="<<rays.size()<<
+	  " pvs="<<mPvsSize<<" pEntropy="<<mPvsEntropy<<" lEntropy="<<mRayLengthEntropy<<endl;
+  };
+  
+  void AddRay(const RayInfo &data) {
+	mValidPvs = false;
+    rays.push_back(data);
+    data.mRay->Ref();
+	if (data.GetRayClass() == RayInfo::PASSING_RAY)
+	  mPassingRays++;
+  }
+	
+  int GetPvsSize() const {
+	return mPvsSize;
+  }
+
+  void SetPvsSize(const int s) {
+	mPvsSize = s;
+	mValidPvs = true;
+  }
+
+  void
+  UpdatePvsSize();
+
+  float
+  ComputePvsEntropy();
+  
+  float
+  ComputeRayLengthEntropy();
+  
+  float
+  ComputeRayTerminationEntropy();
+  
+  void
+  ComputeEntropyImportance();
+
+  void Mail() { mailbox = mailID; }
+  static void NewMail() { mailID++; }
+  bool Mailed() const { return mailbox == mailID; }
+  
+  bool Mailed(const int mail) {
+    return mailbox >= mailID + mail;
+  }
+
+  float GetAvgRayContribution() const {
+	return GetPvsSize()/((float)rays.size() + Limits::Small);
+  }
+
+  float GetImportance() const;
+  
+  float GetSqrRayContribution() const {
+	return sqr(GetPvsSize()/((float)rays.size() + Limits::Small));
+  }
+  
+  // comparator for the 
+  struct less_contribution : public
+  binary_function<const RssTreeLeaf *, const RssTreeLeaf *, bool> {
+	
+	bool operator()(const RssTreeLeaf * a, const RssTreeLeaf *b) {
+	  return a->GetAvgRayContribution() < b->GetAvgRayContribution();
+	}
+  };
+  
+  struct greater_contribution : public
+  binary_function<const RssTreeLeaf *, const RssTreeLeaf *, bool> {
+	
+	bool operator()(const RssTreeLeaf * a, const RssTreeLeaf *b) {
+	  return a->GetAvgRayContribution() > b->GetAvgRayContribution();
+	}
+  };
+  
+  friend bool GreaterContribution(const RssTreeLeaf * a, const RssTreeLeaf *b) {
+	return a->GetAvgRayContribution() > b->GetAvgRayContribution();
+  }
+  
+};
+
+// Inline functions
+inline
+RssTreeNode::RssTreeNode(RssTreeInterior *p):
+  parent(p), axis(-1), depth(p ? p->depth + 1 : 0) {}
+
+
+
+// ---------------------------------------------------------------
+// Main LSDS search class
+// ---------------------------------------------------------------
+class RssTree
+{
+  struct TraversalData
+  {  
+    RssTreeNode *node;
+    AxisAlignedBox3 bbox;
+    int depth;
+    float priority;
+    
+    TraversalData() {}
+
+    TraversalData(RssTreeNode *n, const float p):
+      node(n), priority(p)
+    {}
+
+    TraversalData(RssTreeNode *n,
+				  const AxisAlignedBox3 &b,
+				  const int d):
+      node(n), bbox(b), depth(d) {}
+    
+		
+    // comparator for the 
+    struct less_priority : public
+    binary_function<const TraversalData, const TraversalData, bool> {
+			
+      bool operator()(const TraversalData a, const TraversalData b) {
+		return a.priority < b.priority;
+      }
+      
+    };
+
+    //    ~TraversalData() {}
+    //    TraversalData(const TraversalData &s):node(s.node), bbox(s.bbox), depth(s.depth) {}
+    
+    friend bool operator<(const TraversalData &a,
+						  const TraversalData &b) {
+      //      return a.node->queries.size() < b.node->queries.size();
+      RssTreeLeaf *leafa = (RssTreeLeaf *) a.node;
+      RssTreeLeaf *leafb = (RssTreeLeaf *) b.node;
+#if 0
+	  return
+		leafa->rays.size()*a.bbox.GetVolume()
+		<
+		leafb->rays.size()*b.bbox.GetVolume();
+#endif
+#if 0
+	  return
+		leafa->GetPvsSize()*a.bbox.GetVolume()
+		<
+		leafb->GetPvsSize()*b.bbox.GetVolume();
+#endif
+#if 0
+	  return
+		leafa->GetPvsSize()
+		<
+		leafb->GetPvsSize();
+#endif
+#if 0
+	  return
+		leafa->GetPvsSize()/(leafa->rays.size()+1)
+		>
+		leafb->GetPvsSize()/(leafb->rays.size()+1);
+#endif
+#if 1
+	  return
+		leafa->GetPvsSize()*leafa->rays.size()
+		<
+		leafb->GetPvsSize()*leafb->rays.size();
+#endif
+    }
+  };
+  
+  // simplified data for ray traversal only...
+
+  struct RayTraversalData {
+    
+    RssTreeNode::RayInfo rayData;
+    RssTreeNode *node;
+    
+    RayTraversalData() {}
+    RayTraversalData(RssTreeNode *n,
+					 const RssTreeNode::RayInfo &data):
+      rayData(data), node(n) {}
+  };
+	
+public:
+  /////////////////////////////
+  // The core pointer
+  RssTreeNode *root;
+  
+  /////////////////////////////
+  // Basic properties
+
+  // total number of nodes of the tree
+  int nodes;
+  // axis aligned bounding box of the scene
+  AxisAlignedBox3 bbox;
+
+  // axis aligned bounding box of directions
+  AxisAlignedBox3 dirBBox;
+  
+  /////////////////////////////
+  // Construction parameters
+
+  // epsilon used for the construction
+  float epsilon;
+
+  // ratio between traversal and intersection costs
+  float ct_div_ci;
+  // max depth of the tree
+  int termMaxDepth;
+  // minimal ratio of the volume of the cell and the query volume
+  float termMinSize;
+
+  // minimal pvs per node to still get subdivided
+  int termMinPvs;
+
+  // minimal ray number per node to still get subdivided
+  int termMinRays;
+	
+  // maximal cost ration to subdivide a node
+  float termMaxCostRatio;
+	
+  // maximal contribution per ray to subdivide the node
+  float termMaxRayContribution;
+
+	
+  // randomized construction
+  bool randomize;
+  
+  // type of the splitting to use fo rthe tree construction
+  enum {ESplitRegular, ESplitHeuristic, ESplitHybrid };
+  int splitType;
+
+  bool mSplitUseOnlyDrivingAxis;
+  
+  // use ray space subdivision instead of view space subdivision
+  bool mUseRss;
+
+  // interleave directional and spatial splits based on their costs
+  // if false directional splits are only performed after spatial splits
+  bool mInterleaveDirSplits;
+
+  // depth at which directional splits are performed if mInterleaveDirSplits is false
+  int mDirSplitDepth;
+  
+  // maximal size of the box on which the refdir splitting can be performed
+  // (relative to the scene bbox
+  float refDirBoxMaxSize;
+  
+  // maximum alovable memory in MB
+  float maxTotalMemory;
+
+  // maximum alovable memory for static kd tree in MB
+  float maxStaticMemory;
+
+  // this is used during the construction depending
+  // on the type of the tree and queries...
+  float maxMemory;
+
+
+  // minimal acess time for collapse
+  int accessTimeThreshold;
+
+  // minimal depth at which to perform collapse
+  int minCollapseDepth;
+
+  
+  // reusable array of split candidates
+  vector<SortableEntry> *splitCandidates;
+  /////////////////////////////
+
+  RssStatistics stat;
+	
+  
+  RssTree();
+  virtual ~RssTree();
+
+  virtual void
+  Construct(
+			VssRayContainer &rays,
+			AxisAlignedBox3 *forcedBoundingBox = NULL
+			);
+	
+  // incemental construction
+  virtual void UpdateRays(VssRayContainer &remove,
+						  VssRayContainer &add
+						  );
+
+  virtual void AddRays(
+					   VssRayContainer &add
+					   )
+  {
+	VssRayContainer remove;
+	UpdateRays(remove, add);
+  }
+
+  
+	
+  RssTreeNode *
+  Locate(const Vector3 &v);
+	
+  RssTreeNode *
+  SubdivideNode(RssTreeLeaf *leaf,
+				const AxisAlignedBox3 &box,
+				AxisAlignedBox3 &backBox,
+				AxisAlignedBox3 &frontBox
+				);
+	
+  RssTreeNode *
+  Subdivide(const TraversalData &tdata);
+	
+  int
+  SelectPlane(RssTreeLeaf *leaf,
+			  const AxisAlignedBox3 &box,
+			  float &position,
+			  int &raysBack,
+			  int &raysFront,
+			  int &pvsBack,
+			  int &pvsFront
+			  );
+
+  void
+  SortSplitCandidates(
+					  RssTreeLeaf *node,
+					  const int axis
+					  );
+	
+  
+  // return memory usage in MB
+  float GetMemUsage() const {
+    return 
+      (sizeof(RssTree) +
+       stat.Leaves()*sizeof(RssTreeLeaf) +
+       stat.Interior()*sizeof(RssTreeInterior) +
+       stat.rayRefs*sizeof(RssTreeNode::RayInfo))/(1024.0f*1024.0f);
+  }
+	
+  float GetRayMemUsage() const {
+    return 
+      stat.rays*(sizeof(VssRay))/(1024.0f*1024.0f);
+  }
+  
+
+  float
+  BestCostRatio(
+				RssTreeLeaf *node,
+				int &axis,
+				float &position,
+				int &raysBack,
+				int &raysFront,
+				int &pvsBack,
+				int &pvsFront
+				);
+	
+  float
+  EvalCostRatio(
+				RssTreeLeaf *node,
+				const int axis,
+				const float position,
+				int &raysBack,
+				int &raysFront,
+				int &pvsBack,
+				int &pvsFront
+				);
+
+  float
+  EvalCostRatioHeuristic(
+						 RssTreeLeaf *node,
+						 const int axis,
+						 float &position,
+						 int &raysBack,
+						 int &raysFront,
+						 int &pvsBack,
+						 int &pvsFront
+						 );
+
+  float
+  GetCostRatio(
+			   RssTreeLeaf *leaf,
+			   const int axis,
+			   const float position,
+			   const int raysBack,
+			   const int raysFront,
+			   const int pvsBack,
+			   const int pvsFront
+			   );
+
+  AxisAlignedBox3 GetBBox(const RssTreeNode *node) const {
+    if (node->parent == NULL)
+      return bbox;
+
+    if (!node->IsLeaf())
+      return ((RssTreeInterior *)node)->bbox;
+
+    if (node->parent->axis >= 3)
+      return node->parent->bbox;
+      
+    AxisAlignedBox3 box(node->parent->bbox);
+    if (node->parent->front == node)
+      box.SetMin(node->parent->axis, node->parent->position);
+    else
+      box.SetMax(node->parent->axis, node->parent->position);
+    return box;
+  }
+
+  AxisAlignedBox3 GetShrankedBBox(const RssTreeNode *node) const;
+
+  AxisAlignedBox3 GetDirBBox(const RssTreeNode *node) const {
+
+    if (node->parent == NULL)
+      return dirBBox;
+    
+    if (!node->IsLeaf() )
+      return ((RssTreeInterior *)node)->dirBBox;
+
+    if (node->parent->axis < 3)
+      return node->parent->dirBBox;
+    
+    AxisAlignedBox3 dBBox(node->parent->dirBBox);
+
+    if (node->parent->front == node)
+      dBBox.SetMin(node->parent->axis - 3, node->parent->position);
+    else
+      dBBox.SetMax(node->parent->axis - 3, node->parent->position);
+    return dBBox;
+  }
+  
+  int
+  ReleaseMemory(const int time);
+
+  int
+  CollapseSubtree(RssTreeNode *node, const int time);
+
+  void
+  CountAccess(RssTreeInterior *node, const long time) {
+    node->accesses++;
+    node->lastAccessTime = time;
+  }
+
+  RssTreeNode *
+  SubdivideLeaf(
+				RssTreeLeaf *leaf
+				);
+
+  void
+  RemoveRay(VssRay *ray,
+			vector<RssTreeLeaf *> *affectedLeaves,
+			const bool removeAllScheduledRays
+			);
+
+  //  void
+  //  AddRay(VssRay *ray);
+  void
+  AddRay(RssTreeNode::RayInfo &info);
+  
+  void
+  TraverseInternalNode(
+					   RayTraversalData &data,
+					   stack<RayTraversalData> &tstack);
+
+  void
+  EvaluateLeafStats(const TraversalData &data);
+
+
+  int
+  GetRootPvsSize() const {
+	return GetPvsSize(bbox);
+  }
+  
+  int
+  GetPvsSize(const AxisAlignedBox3 &box) const;
+
+  int
+  CollectPvs(const AxisAlignedBox3 &box,
+			 ObjectContainer &pvs
+			 ) const;
+
+  int
+  CollectRootPvs(
+				 ObjectContainer &pvs
+				 ) const
+  {
+	return CollectPvs(bbox, pvs);
+  }
+  
+	
+  void
+  GetTreeStatistics(
+					float &avgPvsSize,
+					float &avgRayContribution,
+					float &avgPvsEntropy,
+					float &avgRayLengthEntropy,
+					float &avgImportance
+					);
+
+
+  int
+  GenerateRays(const float ratioPerLeaf,
+			   SimpleRayContainer &rays);
+
+  int
+  GenerateRays(const int numberOfRays,
+			   const int numberOfLeaves,
+			   SimpleRayContainer &rays);
+		
+  float
+  GetAvgPvsSize();
+
+  int
+  UpdateSubdivision();
+
+  bool
+  TerminationCriteriaSatisfied(RssTreeLeaf *leaf);
+
+  void
+  CollectLeaves(vector<RssTreeLeaf *> &leaves);
+
+  bool
+  ClipRay(
+		  RssTreeNode::RayInfo &rayInfo,
+		  const AxisAlignedBox3 &box
+		  );
+
+  RssTreeNode *GetRoot() const { return root; }
+
+  bool
+  ValidLeaf(RssTreeLeaf *leaf) const;
+
+  void
+  GenerateLeafRays(RssTreeLeaf *leaf,
+				   const int numberOfRays,
+				   SimpleRayContainer &rays);
+
+
+};
+
+
+#endif // __LSDS_KDTREE_H__
+
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.cpp	(revision 447)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.cpp	(revision 447)
@@ -0,0 +1,2024 @@
+
+// ================================================================
+// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
+// ****************************************************************
+/**
+   The KD tree based LSDS
+*/
+// Initial coding by
+/**
+   @author Jiri Bittner
+*/
+
+// Standard headers
+#include <stack>
+#include <queue>
+#include <algorithm>
+#include <fstream>
+#include <string>
+
+#include "RssTree.h"
+
+#include "Environment.h"
+#include "VssRay.h"
+#include "Intersectable.h"
+#include "Ray.h"
+#include "Containers.h"
+
+#define DEBUG_SPLIT_COST 0
+#define DEBUG_SPLITS 0
+
+// Static variables
+int
+RssTreeLeaf::mailID = 0;
+
+inline void
+AddObject2Pvs(Intersectable *object,
+			  const int side,
+			  int &pvsBack,
+			  int &pvsFront)
+{
+  
+  if (!object)
+	return;
+  
+  if (side <= 0) {
+	if (!object->Mailed() && !object->Mailed(2)) {
+	  pvsBack++;
+	  if (object->Mailed(1))
+		object->Mail(2);
+	  else
+		object->Mail();
+	}
+  }
+  
+  if (side >= 0) {
+	if (!object->Mailed(1) && !object->Mailed(2)) {
+	  pvsFront++;
+	  if (object->Mailed())
+		object->Mail(2);
+	  else
+		object->Mail(1);
+	}
+  }
+}
+
+// Constructor
+RssTree::RssTree()
+{
+  environment->GetIntValue("RssTree.maxDepth", termMaxDepth);
+  environment->GetIntValue("RssTree.minPvs", termMinPvs);
+  environment->GetIntValue("RssTree.minRays", termMinRays);
+  environment->GetFloatValue("RssTree.maxRayContribution", termMaxRayContribution);
+  environment->GetFloatValue("RssTree.maxCostRatio", termMaxCostRatio);
+
+  environment->GetFloatValue("RssTree.minSize", termMinSize);
+  termMinSize = sqr(termMinSize);
+	
+  environment->GetFloatValue("RssTree.refDirBoxMaxSize", refDirBoxMaxSize);
+  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
+  
+  environment->GetFloatValue("RssTree.epsilon", epsilon);
+  environment->GetFloatValue("RssTree.ct_div_ci", ct_div_ci);
+	
+  environment->GetFloatValue("RssTree.maxTotalMemory", maxTotalMemory);
+  environment->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
+  
+  
+
+
+  float refDirAngle;
+  environment->GetFloatValue("RssTree.refDirAngle", refDirAngle);
+  
+  environment->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
+  //= 1000;
+  environment->GetIntValue("RssTree.minCollapseDepth", minCollapseDepth);
+  //  int minCollapseDepth = 4;
+
+  //  pRefDirThresh = cos(0.5*M_PI - M_PI*refDirAngle/180.0);
+  //  cosRefDir = cos(M_PI*refDirAngle/180.0);
+  //  sinRefDir = sin(M_PI*refDirAngle/180.0);
+  
+  
+  // split type
+  char sname[128];
+  environment->GetStringValue("RssTree.splitType", sname);
+  string name(sname);
+	
+  if (name.compare("regular") == 0)
+    splitType = ESplitRegular;
+  else
+    if (name.compare("heuristic") == 0)
+      splitType = ESplitHeuristic;
+	else
+	  if (name.compare("hybrid") == 0)
+		splitType = ESplitHybrid;
+	  else {
+		cerr<<"Invalid RssTree split type "<<name<<endl;
+		exit(1);
+	  }
+	
+  environment->GetBoolValue("RssTree.randomize", randomize);
+  environment->GetBoolValue("RssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
+  environment->GetBoolValue("RssTree.useRss", mUseRss);
+
+  environment->GetBoolValue("RssTree.interleaveDirSplits", mInterleaveDirSplits);
+  environment->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);
+
+  root = NULL;
+  
+  splitCandidates = new vector<SortableEntry>;
+}
+
+
+RssTree::~RssTree()
+{
+  if (root)
+    delete root;
+}
+
+
+
+
+void
+RssStatistics::Print(ostream &app) const
+{
+  app << "===== RssTree statistics ===============\n";
+
+  app << "#N_RAYS ( Number of rays )\n"
+      << rays <<endl;
+
+  app << "#N_INITPVS ( Initial PVS size )\n"
+      << initialPvsSize <<endl;
+  
+  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
+
+  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
+
+  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
+  for (int i=0; i<7; i++)
+    app << splits[i] <<" ";
+  app <<endl;
+
+  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
+    rayRefs << "\n";
+
+  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
+    rayRefs/(double)rays << "\n";
+
+  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
+    rayRefs/(double)Leaves() << "\n";
+
+  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
+    maxRayRefs << "\n";
+
+
+  //  app << setprecision(4);
+
+  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
+    maxDepthNodes*100/(double)Leaves()<<endl;
+
+  app << "#N_PMINPVSLEAVES  ( Percentage of leaves with minCost )\n"<<
+    minPvsNodes*100/(double)Leaves()<<endl;
+
+  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\n"<<
+    minRaysNodes*100/(double)Leaves()<<endl;
+	
+  app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
+    minSizeNodes*100/(double)Leaves()<<endl;
+
+  app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
+    maxRayContribNodes*100/(double)Leaves()<<endl;
+
+  app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
+    maxCostRatioNodes*100/(double)Leaves()<<endl;
+
+  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
+    addedRayRefs<<endl;
+
+  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
+    removedRayRefs<<endl;
+
+  //  app << setprecision(4);
+
+  app << "#N_CTIME  ( Construction time [s] )\n"
+      << Time() << " \n";
+
+  app << "===== END OF RssTree statistics ==========\n";
+
+}
+
+
+void
+RssTreeLeaf::UpdatePvsSize()
+{
+  if (!mValidPvs) {
+	Intersectable::NewMail();
+	int pvsSize = 0;
+	for(RssTreeNode::RayInfoContainer::iterator ri = rays.begin();
+		ri != rays.end();
+		ri++)
+	  if ((*ri).mRay->IsActive()) {
+		Intersectable *object;
+#if BIDIRECTIONAL_RAY
+		object = (*ri).mRay->mOriginObject;
+		if (object && !object->Mailed()) {
+		  pvsSize++;
+		  object->Mail();
+		}
+#endif
+		object = (*ri).mRay->mTerminationObject;
+		if (object && !object->Mailed()) {
+		  pvsSize++;
+		  object->Mail();
+		}
+	  }
+	ComputeEntropyImportance();
+
+	mPvsSize = pvsSize;
+	mValidPvs = true;
+  }
+}
+
+bool
+RssTree::ClipRay(
+				 RssTreeNode::RayInfo &rayInfo,
+				 const AxisAlignedBox3 &box
+				 )
+{
+  float tmin, tmax;
+  static Ray ray;
+  ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
+  box.ComputeMinMaxT(ray, &tmin, &tmax);
+  if (tmin >= tmax)
+	return false;
+  
+  // now check if the ray origin lies inside the box
+  if ( tmax < rayInfo.mRay->GetSize() ) {
+	// this ray does not leave the box
+	rayInfo.mRay->SetupEndPoints(
+								 ray.Extrap(tmax),
+								 rayInfo.mRay->mTermination
+								 );
+	return true;
+  }
+
+  return false;
+}
+
+
+void
+RssTree::Construct(
+				   VssRayContainer &rays,
+				   // forced bounding box is only used when computing from-box
+				   // visibility
+				   AxisAlignedBox3 *forcedBoundingBox
+				   )
+{
+  stat.Start();
+  
+  maxMemory = maxStaticMemory;
+
+  if (root)
+    delete root;
+	
+  root = new RssTreeLeaf(NULL, rays.size());
+  // first construct a leaf that will get subdivide
+  RssTreeLeaf *leaf = (RssTreeLeaf *) root;
+
+  stat.nodes = 1;
+  
+  bbox.Initialize();
+  dirBBox.Initialize();
+
+  if (mUseRss)
+	forcedBoundingBox = NULL;
+	
+  for(VssRayContainer::const_iterator ri = rays.begin();
+      ri != rays.end();
+      ri++) {
+		
+	RssTreeNode::RayInfo info(*ri);
+	if (forcedBoundingBox)
+	  if (!ClipRay(info, *forcedBoundingBox))
+		continue;
+
+	leaf->AddRay(info);
+		
+    bbox.Include((*ri)->GetOrigin());
+    bbox.Include((*ri)->GetTermination());
+    
+		
+	dirBBox.Include(Vector3(
+							(*ri)->GetDirParametrization(0),
+							(*ri)->GetDirParametrization(1),
+							0
+							)
+					);
+  }
+	
+	
+  if ( forcedBoundingBox ) 
+	bbox = *forcedBoundingBox;
+	
+  cout<<"Bbox = "<<bbox<<endl;
+  cout<<"Dirr Bbox = "<<dirBBox<<endl;
+
+  stat.rays = leaf->rays.size();
+  leaf->UpdatePvsSize();
+  leaf->ComputeEntropyImportance();
+
+  stat.initialPvsSize = leaf->GetPvsSize();
+  // Subdivide();
+  root = Subdivide(TraversalData(leaf, bbox, 0));
+
+  if (splitCandidates) {
+    // force realease of this vector
+    delete splitCandidates;
+    splitCandidates = new vector<SortableEntry>;
+  }
+  
+  stat.Stop();
+
+  stat.Print(cout);
+  cout<<"#Total memory="<<GetMemUsage()<<endl;
+
+}
+
+int
+RssTree::UpdateSubdivision()
+{
+  priority_queue<TraversalData> tStack;
+  //  stack<TraversalData> tStack;
+  
+  tStack.push(TraversalData(root, bbox, 0));
+	
+  AxisAlignedBox3 backBox;
+  AxisAlignedBox3 frontBox;
+
+  maxMemory = maxTotalMemory;
+  int subdivided = 0;
+  int lastMem = 0;
+  while (!tStack.empty()) {
+		
+	float mem = GetMemUsage();
+		
+	if ( lastMem/10 != ((int)mem)/10) {
+	  cout<<mem<<" MB"<<endl;
+	}
+	lastMem = (int)mem;
+		
+	if (  mem > maxMemory ) {
+      // count statistics on unprocessed leafs
+      while (!tStack.empty()) {
+		//				EvaluateLeafStats(tStack.top());
+		tStack.pop();
+      }
+      break;
+    }
+    
+    TraversalData data = tStack.top();
+    tStack.pop();
+
+	if (data.node->IsLeaf()) {
+	  RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
+										data.bbox,
+										backBox,
+										frontBox
+										);
+	  if (!node->IsLeaf()) {
+		subdivided++;
+		RssTreeInterior *interior = (RssTreeInterior *) node;
+		// push the children on the stack
+		tStack.push(TraversalData(interior->back, backBox, data.depth+1));
+		tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
+	  } else {
+		//      EvaluateLeafStats(data);
+	  }
+	} else {
+	  RssTreeInterior *interior = (RssTreeInterior *) data.node;
+	  tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
+	  tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
+	}
+  }
+  return subdivided;
+}
+
+
+RssTreeNode *
+RssTree::Subdivide(const TraversalData &tdata)
+{
+  RssTreeNode *result = NULL;
+
+  priority_queue<TraversalData> tStack;
+  //  stack<TraversalData> tStack;
+  
+  tStack.push(tdata);
+
+  AxisAlignedBox3 backBox;
+  AxisAlignedBox3 frontBox;
+
+  
+  int lastMem = 0;
+  while (!tStack.empty()) {
+
+	float mem = GetMemUsage();
+		
+	if ( lastMem/10 != ((int)mem)/10) {
+	  cout<<mem<<" MB"<<endl;
+	}
+	lastMem = (int)mem;
+		
+	if (  mem > maxMemory ) {
+      // count statistics on unprocessed leafs
+      while (!tStack.empty()) {
+		EvaluateLeafStats(tStack.top());
+		tStack.pop();
+      }
+      break;
+    }
+    
+    TraversalData data = tStack.top();
+    tStack.pop();
+
+#if DEBUG_SPLITS
+	Debug<<"#Splitting node"<<endl;
+	data.node->Print(Debug);
+#endif
+	RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
+									  data.bbox,
+									  backBox,
+									  frontBox
+									  );
+	
+
+	if (result == NULL)
+      result = node;
+    
+    if (!node->IsLeaf()) {
+			
+      RssTreeInterior *interior = (RssTreeInterior *) node;
+      // push the children on the stack
+      tStack.push(TraversalData(interior->back, backBox, data.depth+1));
+      tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
+
+#if DEBUG_SPLITS
+	  Debug<<"#New nodes"<<endl;
+	  interior->back->Print(Debug);
+	  interior->front->Print(Debug);
+	  Debug<<"#####################################"<<endl;
+#endif
+	  
+    } else {
+      EvaluateLeafStats(data);
+    }
+  }
+
+  return result;
+}
+
+
+// returns selected plane for subdivision
+int
+RssTree::SelectPlane(
+					 RssTreeLeaf *leaf,
+					 const AxisAlignedBox3 &box,
+					 float &position,
+					 int &raysBack,
+					 int &raysFront,
+					 int &pvsBack,
+					 int &pvsFront
+					 )
+{
+
+  int minDirDepth = 6;
+  int axis;
+  float costRatio;
+	
+  costRatio = BestCostRatio(leaf,
+							axis,
+							position,
+							raysBack,
+							raysFront,
+							pvsBack,
+							pvsFront
+							);
+#if DEBUG_SPLIT_COST
+  cout<<axis<<" r="<<costRatio<<endl;
+#endif
+  if (costRatio > termMaxCostRatio) {
+	//		cout<<"Too big cost ratio "<<costRatio<<endl;
+	stat.maxCostRatioNodes++;
+	return -1;
+  }
+	
+#if 0	
+  cout<<
+	"pvs="<<leaf->mPvsSize<<
+	" rays="<<leaf->rays.size()<<
+	" rc="<<leaf->GetAvgRayContribution()<<
+	" axis="<<axis<<endl;
+#endif
+	
+  return axis;
+}
+
+
+float
+RssTree::GetCostRatio(
+					  RssTreeLeaf *leaf,
+					  const int axis,
+					  const float position,
+					  const int raysBack,
+					  const int raysFront,
+					  const int pvsBack,
+					  const int pvsFront
+					  )
+{
+  bool costImportance = true;
+
+  float ratio;
+  AxisAlignedBox3 box;
+  float minBox, maxBox;
+
+  if (axis < 3) {
+	box = GetBBox(leaf);
+	minBox = box.Min(axis);
+	maxBox = box.Max(axis);
+  }	else {
+	box = GetDirBBox(leaf);
+	minBox = box.Min(axis-3);
+	maxBox = box.Max(axis-3);
+  }
+	
+  float sizeBox = maxBox - minBox;
+
+  int pvsSize = leaf->GetPvsSize();
+
+  if (!costImportance) {
+	//		float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+	float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
+	float newCost = ct_div_ci + sum/sizeBox;
+	float oldCost = pvsSize;
+	ratio = newCost/oldCost;
+  } else {
+	// importance based cost
+#if 0
+	float newContrib =
+	  ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
+	   (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
+		
+	//			float newContrib =
+	//				sqr(pvsBack/(raysBack + Limits::Small)) +
+	//				sqr(pvsFront/(raysFront + Limits::Small));
+	float oldContrib = sqr(leaf->GetAvgRayContribution());
+	ratio = oldContrib/newContrib;
+#else
+#if 1
+	float newCost = raysBack*pvsBack  + raysFront*pvsFront;
+	float oldCost = leaf->rays.size()*pvsSize;
+	ratio = newCost/oldCost;
+#else
+#if 0
+	float newCost = (pvsBack  + pvsFront)*0.5f;
+	float oldCost = pvsSize;
+	ratio = newCost/oldCost;
+#else
+	float newCost = abs(raysBack  - raysFront);
+	float oldCost = leaf->rays.size();
+	ratio = newCost/oldCost;
+#endif
+#endif
+#endif
+  }
+	
+  return ratio;
+}
+							
+
+float
+RssTree::EvalCostRatio(
+					   RssTreeLeaf *leaf,
+					   const int axis,
+					   const float position,
+					   int &raysBack,
+					   int &raysFront,
+					   int &pvsBack,
+					   int &pvsFront
+					   )
+{
+  raysBack = 0;
+  raysFront = 0;
+  pvsFront = 0;
+  pvsBack = 0;
+
+	
+  Intersectable::NewMail(3);
+	
+  if (axis <= RssTreeNode::SPLIT_Z) {
+    // this is the main ray classification loop!
+    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		ri != leaf->rays.end();
+		ri++)
+      if ((*ri).mRay->IsActive()) {
+				
+		// determine the side of this ray with respect to the plane
+		int side = (*ri).ComputeRaySide(axis, position);
+		//				(*ri).mRay->mSide = side;
+				
+		if (side <= 0)
+		  raysBack++;
+				
+		if (side >= 0)
+		  raysFront++;
+				
+		AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+      }
+		
+  } else {
+		
+	// directional split
+    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		ri != leaf->rays.end();
+		ri++)
+      if ((*ri).mRay->IsActive()) {
+				
+		// determine the side of this ray with respect to the plane
+		int side;
+		if ((*ri).mRay->GetDirParametrization(axis - 3) <= position)
+		  side = -1;
+		else
+		  side = 1;
+		
+		if (side <= 0)
+		  raysBack++;
+				
+		if (side >= 0)
+		  raysFront++;
+
+		//				(*ri).mRay->mSide = side;
+		AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+				
+      }
+  }
+
+  float ratio = GetCostRatio(
+							 leaf,
+							 axis,
+							 position,
+							 raysBack,
+							 raysFront,
+							 pvsBack,
+							 pvsFront);
+
+  //	cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
+  //  float oldCost = leaf->rays.size();
+
+  //	cout<<"ratio="<<ratio<<endl;
+
+  return ratio;
+}
+
+float
+RssTree::BestCostRatio(
+					   RssTreeLeaf *leaf,
+					   int &axis,
+					   float &position,
+					   int &raysBack,
+					   int &raysFront,
+					   int &pvsBack,
+					   int &pvsFront
+					   )
+{
+  int nRaysBack[6], nRaysFront[6];
+  int nPvsBack[6], nPvsFront[6];
+  float nPosition[6];
+  float nCostRatio[6];
+  int bestAxis = -1;
+	
+  AxisAlignedBox3 sBox = GetBBox(leaf);
+  AxisAlignedBox3 dBox = GetDirBBox(leaf);
+  // int sAxis = box.Size().DrivingAxis();
+  int sAxis = sBox.Size().DrivingAxis();
+  int dAxis = dBox.Size().DrivingAxis() + 3;
+
+
+  float dirSplitBoxSize = 0.01f;
+  bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
+		
+  
+  for (axis = 0; axis < 5; axis++)
+	if (
+		(axis < 3 && (leaf->depth < mDirSplitDepth ||  mInterleaveDirSplits)) ||
+		(axis >= 3 && (leaf->depth >= mDirSplitDepth))
+		) {
+	  if (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis) {
+		
+		if (splitType == ESplitRegular) {
+		  if (axis < 3)
+			nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
+		  else
+			nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
+		  
+		  nCostRatio[axis] = EvalCostRatio(leaf,
+										   axis,
+										   nPosition[axis],
+										   nRaysBack[axis],
+										   nRaysFront[axis],
+										   nPvsBack[axis],
+										   nPvsFront[axis]
+										   );
+		} else
+		  if (splitType == ESplitHeuristic) {
+			nCostRatio[axis] = EvalCostRatioHeuristic(
+													  leaf,
+													  axis,
+													  nPosition[axis],
+													  nRaysBack[axis],
+													  nRaysFront[axis],
+													  nPvsBack[axis],
+													  nPvsFront[axis]);
+		  } else
+			if (splitType == ESplitHybrid) {
+			  if (leaf->depth > 7)
+				nCostRatio[axis] = EvalCostRatioHeuristic(
+														  leaf,
+														  axis,
+														  nPosition[axis],
+														  nRaysBack[axis],
+														  nRaysFront[axis],
+														  nPvsBack[axis],
+														  nPvsFront[axis]);
+			  else {
+				if (axis < 3)
+				  nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
+				else
+				  nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
+				
+				nCostRatio[axis] = EvalCostRatio(leaf,
+												 axis,
+												 nPosition[axis],
+												 nRaysBack[axis],
+												 nRaysFront[axis],
+												 nPvsBack[axis],
+												 nPvsFront[axis]
+												 );
+			  }
+			} else {
+			  cerr<<"RssTree: Unknown split heuristics\n";
+			  exit(1);
+			}
+		
+		
+		if ( bestAxis == -1)
+		  bestAxis = axis;
+		else
+		  if ( nCostRatio[axis] < nCostRatio[bestAxis] )
+			bestAxis = axis;
+	  }
+	}
+  
+  axis = bestAxis;
+  position = nPosition[bestAxis];
+
+  raysBack = nRaysBack[bestAxis];
+  raysFront = nRaysFront[bestAxis];
+
+  pvsBack = nPvsBack[bestAxis];
+  pvsFront = nPvsFront[bestAxis];
+	
+  return nCostRatio[bestAxis];
+}
+
+	
+float
+RssTree::EvalCostRatioHeuristic(
+								RssTreeLeaf *leaf,
+								const int axis,
+								float &bestPosition,
+								int &raysBack,
+								int &raysFront,
+								int &pvsBack,
+								int &pvsFront
+								)
+{
+  AxisAlignedBox3 box;
+  float minBox, maxBox;
+	
+  if (axis < 3) {
+	box = GetBBox(leaf);
+	minBox = box.Min(axis);
+	maxBox = box.Max(axis);
+  } else {
+	box = GetDirBBox(leaf);
+	minBox = box.Min(axis-3);
+	maxBox = box.Max(axis-3);
+  }
+	
+  SortSplitCandidates(leaf, axis);
+  
+  // go through the lists, count the number of objects left and right
+  // and evaluate the following cost funcion:
+  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
+	
+  int rl=0, rr = leaf->rays.size();
+  int pl=0, pr = leaf->GetPvsSize();
+  float sizeBox = maxBox - minBox;
+	
+  float minBand = minBox + 0.1*(maxBox - minBox);
+  float maxBand = minBox + 0.9*(maxBox - minBox);
+	
+  float minRatio = 1e20;
+	
+  Intersectable::NewMail();
+  // set all object as belonging to the fron pvs
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+	  ri != leaf->rays.end();
+	  ri++)
+	if ((*ri).mRay->IsActive()) {
+	  Intersectable *object = (*ri).mRay->mTerminationObject;
+	  if (object)
+		if (!object->Mailed()) {
+		  object->Mail();
+		  object->mCounter = 1;
+		} else
+		  object->mCounter++;
+	}
+	
+  Intersectable::NewMail();
+	
+  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
+	  ci < splitCandidates->end();
+	  ci++) {
+	VssRay *ray;
+	switch ((*ci).type) {
+	case SortableEntry::ERayMin: {
+	  rr--;
+	  rl++;
+	  ray = (VssRay *) (*ci).data;
+	  Intersectable *object = ray->mTerminationObject;
+	  if (object) {
+		if (!object->Mailed()) {
+		  object->Mail();
+		  pl++;
+		}
+		if (--object->mCounter == 0)
+		  pr--;
+	  }
+	  break;
+	}
+	}
+
+	float position = (*ci).value;
+	
+	if (position > minBand && position < maxBand) {
+			
+	  float ratio = GetCostRatio(
+								 leaf,
+								 axis,
+								 position,
+								 rl,
+								 rr,
+								 pl,
+								 pr);
+			
+			
+	  //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
+	  //      cout<<"cost= "<<sum<<endl;
+			
+	  if (ratio < minRatio) {
+		minRatio = ratio;
+		bestPosition = position;
+				
+		raysBack = rl;
+		raysFront = rr;
+				
+		pvsBack = pl;
+		pvsFront = pr;
+				
+      }
+    }
+  }
+
+  
+  //  cout<<"===================="<<endl;
+  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
+  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
+  return minRatio;
+}
+
+void
+RssTree::SortSplitCandidates(
+							 RssTreeLeaf *node,
+							 const int axis
+							 )
+{
+  
+  splitCandidates->clear();
+  
+  int requestedSize = 2*(node->rays.size());
+  // creates a sorted split candidates array
+  if (splitCandidates->capacity() > 500000 &&
+      requestedSize < (int)(splitCandidates->capacity()/10) ) {
+    
+    delete splitCandidates;
+    splitCandidates = new vector<SortableEntry>;
+  }
+  
+  splitCandidates->reserve(requestedSize);
+
+  // insert all queries 
+  for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
+      ri < node->rays.end();
+      ri++) {
+	if ((*ri).mRay->IsActive()) {
+	  if (axis < 3) {
+		splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
+												 (*ri).ExtrapOrigin(axis),
+												 (void *)(*ri).mRay)
+								   );
+	  } else {
+		float pos = (*ri).mRay->GetDirParametrization(axis-3);
+		splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
+												 pos,
+												 (void *)(*ri).mRay)
+								   );
+	  }
+	}
+  }
+  
+  stable_sort(splitCandidates->begin(), splitCandidates->end());
+}
+
+
+void
+RssTree::EvaluateLeafStats(const TraversalData &data)
+{
+
+  // the node became a leaf -> evaluate stats for leafs
+  RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
+
+  if (data.depth >= termMaxDepth)
+    stat.maxDepthNodes++;
+  
+  //  if ( (int)(leaf->rays.size()) < termMinCost)
+  //    stat.minCostNodes++;
+  if ( leaf->GetPvsSize() < termMinPvs)
+	stat.minPvsNodes++;
+
+  if ( leaf->GetPvsSize() < termMinRays)
+	stat.minRaysNodes++;
+
+  if (leaf->GetAvgRayContribution() > termMaxRayContribution )
+	stat.maxRayContribNodes++;
+	
+  if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
+	stat.minSizeNodes++;
+  }
+
+  if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
+    stat.maxRayRefs = leaf->rays.size();
+
+}
+
+bool
+RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
+{
+  return ( (leaf->GetPvsSize() < termMinPvs) ||
+		   (leaf->rays.size() < termMinRays) ||
+		   //			 (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
+		   (leaf->depth >= termMaxDepth) ||
+		   (SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize) ||
+		   (mUseRss && leaf->mPassingRays == leaf->rays.size())
+		   );
+}
+
+
+RssTreeNode *
+RssTree::SubdivideNode(
+					   RssTreeLeaf *leaf,
+					   const AxisAlignedBox3 &box,
+					   AxisAlignedBox3 &backBBox,
+					   AxisAlignedBox3 &frontBBox
+					   )
+{
+  
+  if (TerminationCriteriaSatisfied(leaf)) {
+#if 0
+	if (leaf->depth >= termMaxDepth) {
+	  cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
+	  cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
+	}
+#endif
+		
+	return leaf;
+  }
+	
+  float position;
+	
+  // first count ray sides
+  int raysBack;
+  int raysFront;
+  int pvsBack;
+  int pvsFront;
+	
+  // select subdivision axis
+  int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
+  Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
+  
+  if (axis == -1) {
+    return leaf;
+  }
+  
+  stat.nodes+=2;
+  stat.splits[axis]++;
+
+  // add the new nodes to the tree
+  RssTreeInterior *node = new RssTreeInterior(leaf->parent);
+
+  node->axis = axis;
+  node->position = position;
+  node->bbox = box;
+  node->dirBBox = GetDirBBox(leaf);
+  
+  backBBox = box;
+  frontBBox = box;
+
+  RssTreeLeaf *back = new RssTreeLeaf(node, raysBack);
+  RssTreeLeaf *front = new RssTreeLeaf(node, raysFront);
+
+  // replace a link from node's parent
+  if (  leaf->parent )
+    leaf->parent->ReplaceChildLink(leaf, node);
+  // and setup child links
+  node->SetupChildLinks(back, front);
+	
+  if (axis <= RssTreeNode::SPLIT_Z) {
+	backBBox.SetMax(axis, position);
+    frontBBox.SetMin(axis, position);
+		
+	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		ri != leaf->rays.end();
+		ri++) {
+      if ((*ri).mRay->IsActive()) {
+				
+		// first unref ray from the former leaf
+		(*ri).mRay->Unref();
+
+		// Debug << "computed t: " << (*ri).mRay->mT << endl;
+		// determine the side of this ray with respect to the plane
+		int side = node->ComputeRaySide(*ri);
+		
+		if (side == 0) {
+		  if ((*ri).mRay->HasPosDir(axis)) {
+			back->AddRay(*ri);
+			front->AddRay(*ri);
+		  } else {
+			back->AddRay(*ri);
+			front->AddRay(*ri);
+		  }
+		} else
+		  if (side == 1) 
+			front->AddRay(*ri);
+		  else
+			back->AddRay(*ri);
+      } else
+		(*ri).mRay->Unref();
+    }
+  } else {
+    // rays front/back
+    
+    
+    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		ri != leaf->rays.end();
+		ri++) {
+      if ((*ri).mRay->IsActive()) {
+		// first unref ray from the former leaf
+		(*ri).mRay->Unref();
+
+		int side;
+		if ((*ri).mRay->GetDirParametrization(axis - 3) <= position)
+		  side = -1;
+		else
+		  side = 1;
+				
+		if (side == 1)
+		  front->AddRay(*ri);
+		else
+		  back->AddRay(*ri);
+				
+      } else
+		(*ri).mRay->Unref();
+    }
+  }
+
+#if 0
+  front->SetPvsSize(pvsFront);
+  back->SetPvsSize(pvsBack);
+  // compute entropy as well
+  front->ComputeEntropyImportance();
+  back->ComputeEntropyImportance();
+#else
+  front->UpdatePvsSize();
+  back->UpdatePvsSize();
+#endif
+
+  
+  // update stats
+  stat.rayRefs -= leaf->rays.size();
+  stat.rayRefs += raysBack + raysFront;
+
+  
+  delete leaf;
+  return node;
+}
+
+
+
+
+
+
+int
+RssTree::ReleaseMemory(const int time)
+{
+  stack<RssTreeNode *> tstack;
+  
+  // find a node in the tree which subtree will be collapsed
+  int maxAccessTime = time - accessTimeThreshold;
+  int released;
+
+  tstack.push(root);
+
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
+    if (!node->IsLeaf()) {
+      RssTreeInterior *in = (RssTreeInterior *)node;
+      //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
+      if (in->depth >= minCollapseDepth &&
+		  in->lastAccessTime <= maxAccessTime) {
+		released = CollapseSubtree(node, time);
+		break;
+      }
+      
+      if (in->back->GetAccessTime() <  
+		  in->front->GetAccessTime()) {
+		tstack.push(in->front);
+		tstack.push(in->back);
+      } else {
+		tstack.push(in->back);
+		tstack.push(in->front);
+      }
+    }
+  }
+
+  while (tstack.empty()) {
+    // could find node to collaps...
+    //    cout<<"Could not find a node to release "<<endl;
+    break;
+  }
+  
+  return released;
+}
+
+
+
+
+RssTreeNode *
+RssTree::SubdivideLeaf(
+					   RssTreeLeaf *leaf
+					   )
+{
+  RssTreeNode *node = leaf;
+	
+  AxisAlignedBox3 leafBBox = GetBBox(leaf);
+
+  static int pass = 0;
+  pass ++;
+	
+  // check if we should perform a dynamic subdivision of the leaf
+  if (!TerminationCriteriaSatisfied(leaf)) {
+    
+	// memory check and realese...
+    if (GetMemUsage() > maxTotalMemory) {
+      ReleaseMemory( pass );
+    }
+    
+    AxisAlignedBox3 backBBox, frontBBox;
+
+    // subdivide the node
+    node =
+      SubdivideNode(leaf,
+					leafBBox,
+					backBBox,
+					frontBBox
+					);
+  }
+	
+  return node;
+}
+
+
+
+void
+RssTree::UpdateRays(VssRayContainer &remove,
+					VssRayContainer &add
+					)
+{
+  RssTreeLeaf::NewMail();
+
+  // schedule rays for removal
+  for(VssRayContainer::const_iterator ri = remove.begin();
+      ri != remove.end();
+      ri++) {
+    (*ri)->ScheduleForRemoval();
+  }
+
+  int inactive=0;
+
+  for(VssRayContainer::const_iterator ri = remove.begin();
+      ri != remove.end();
+      ri++) {
+    if ((*ri)->ScheduledForRemoval())
+	  //      RemoveRay(*ri, NULL, false);
+	  // !!! BUG - with true it does not work correctly - aggreated delete
+      RemoveRay(*ri, NULL, true);
+    else
+      inactive++;
+  }
+
+
+  //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
+  
+  for(VssRayContainer::const_iterator ri = add.begin();
+      ri != add.end();
+      ri++) {
+	RssTreeNode::RayInfo info(*ri);
+	if (ClipRay(info, bbox))
+	  AddRay(info);
+  }
+}
+
+
+void
+RssTree::RemoveRay(VssRay *ray,
+				   vector<RssTreeLeaf *> *affectedLeaves,
+				   const bool removeAllScheduledRays
+				   )
+{
+	
+  stack<RayTraversalData> tstack;
+
+  tstack.push(RayTraversalData(root, RssTreeNode::RayInfo(ray)));
+  
+  RayTraversalData data;
+
+  // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
+
+  while (!tstack.empty()) {
+    data = tstack.top();
+    tstack.pop();
+
+    if (!data.node->IsLeaf()) {
+      // split the set of rays in two groups intersecting the
+      // two subtrees
+
+      TraverseInternalNode(data, tstack);
+      
+    } else {
+      // remove the ray from the leaf
+      // find the ray in the leaf and swap it with the last ray...
+      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
+      
+      if (!leaf->Mailed()) {
+		leaf->Mail();
+		if (affectedLeaves)
+		  affectedLeaves->push_back(leaf);
+	
+		if (removeAllScheduledRays) {
+		  int tail = leaf->rays.size()-1;
+
+		  for (int i=0; i < (int)(leaf->rays.size()); i++) {
+			if (leaf->rays[i].mRay->ScheduledForRemoval()) {
+			  // find a ray to replace it with
+			  while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
+				stat.removedRayRefs++;
+				leaf->rays[tail].mRay->Unref();
+				leaf->rays.pop_back();
+				tail--;
+			  }
+
+			  if (tail < i)
+				break;
+	      
+			  stat.removedRayRefs++;
+			  leaf->rays[i].mRay->Unref();
+			  leaf->rays[i] = leaf->rays[tail];
+			  leaf->rays.pop_back();
+			  tail--;
+			}
+		  }
+		}
+      }
+      
+      if (!removeAllScheduledRays)
+		for (int i=0; i < (int)leaf->rays.size(); i++) {
+		  if (leaf->rays[i].mRay == ray) {
+			stat.removedRayRefs++;
+			ray->Unref();
+			leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
+			leaf->rays.pop_back();
+			// check this ray again
+			break;
+		  }
+		}
+      
+    }
+  }
+
+  if (ray->RefCount() != 0) {
+    cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
+    exit(1);
+  }
+  
+}
+
+
+void
+RssTree::AddRay(RssTreeNode::RayInfo &info)
+{
+
+  stack<RayTraversalData> tstack;
+  
+  tstack.push(RayTraversalData(root, info));
+  
+  RayTraversalData data;
+  
+  while (!tstack.empty()) {
+    data = tstack.top();
+    tstack.pop();
+
+    if (!data.node->IsLeaf()) {
+      TraverseInternalNode(data, tstack);
+    } else {
+      // remove the ray from the leaf
+      // find the ray in the leaf and swap it with the last ray...
+      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
+      leaf->AddRay(data.rayData);
+      stat.addedRayRefs++;
+    }
+  }
+}
+
+void
+RssTree::TraverseInternalNode(
+							  RayTraversalData &data,
+							  stack<RayTraversalData> &tstack)
+{
+  RssTreeInterior *in =  (RssTreeInterior *) data.node;
+  
+  if (in->axis <= RssTreeNode::SPLIT_Z) {
+
+    // determine the side of this ray with respect to the plane
+    int side = in->ComputeRaySide(data.rayData
+								  );
+    
+  
+    if (side == 0) {
+      if (data.rayData.mRay->HasPosDir(in->axis)) {
+		tstack.push(RayTraversalData(in->back,
+									 data.rayData)
+					);
+				
+		tstack.push(RayTraversalData(in->front,
+									 data.rayData)
+					);
+	
+      } else {
+		tstack.push(RayTraversalData(in->back,
+									 data.rayData
+									 )
+					);
+				
+		tstack.push(RayTraversalData(in->front,
+									 data.rayData)
+					);
+				
+				
+      }
+    } else
+      if (side == 1)
+		tstack.push(RayTraversalData(in->front, data.rayData));
+      else
+		tstack.push(RayTraversalData(in->back, data.rayData));
+  } 
+  else {
+    // directional split
+	if (data.rayData.mRay->GetDirParametrization(in->axis - 3) <= in->position)
+	  tstack.push(RayTraversalData(in->back, data.rayData));
+	else
+	  tstack.push(RayTraversalData(in->front, data.rayData));
+  }
+}
+
+
+int
+RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
+{
+  // first count all rays in the subtree
+  // use mail 1 for this purpose
+  stack<RssTreeNode *> tstack;
+  int rayCount = 0;
+  int totalRayCount = 0;
+  int collapsedNodes = 0;
+
+#if DEBUG_COLLAPSE
+  cout<<"Collapsing subtree"<<endl;
+  cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
+  cout<<"depth="<<(int)sroot->depth<<endl;
+#endif
+
+  //    tstat.collapsedSubtrees++;
+  //    tstat.collapseDepths += (int)sroot->depth;
+  //    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
+  
+  tstack.push(sroot);
+  VssRay::NewMail();
+	
+  while (!tstack.empty()) {
+    collapsedNodes++;
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+
+    if (node->IsLeaf()) {
+      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
+      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		  ri != leaf->rays.end();
+		  ri++) {
+				
+		totalRayCount++;
+		if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
+		  (*ri).mRay->Mail();
+		  rayCount++;
+		}
+      }
+    } else {
+      tstack.push(((RssTreeInterior *)node)->back);
+      tstack.push(((RssTreeInterior *)node)->front);
+    }
+  }
+  
+  VssRay::NewMail();
+
+  // create a new node that will hold the rays
+  RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
+  if (  newLeaf->parent )
+    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
+  
+
+  tstack.push( sroot );
+
+  while (!tstack.empty()) {
+
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+
+    if (node->IsLeaf()) {
+      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
+      
+      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		  ri != leaf->rays.end();
+		  ri++) {
+				
+		// unref this ray from the old node
+				
+		if ((*ri).mRay->IsActive()) {
+		  (*ri).mRay->Unref();
+		  if (!(*ri).mRay->Mailed()) {
+			(*ri).mRay->Mail();
+			newLeaf->AddRay(*ri);
+		  }
+		} else
+		  (*ri).mRay->Unref();
+				
+      }
+    } else {
+      tstack.push(((RssTreeInterior *)node)->back);
+      tstack.push(((RssTreeInterior *)node)->front);
+    }
+  }
+  
+  // delete the node and all its children
+  delete sroot;
+  
+  //   for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
+  //       ri != newLeaf->rays.end();
+  //       ri++)
+  //     (*ri).ray->UnMail(2);
+
+
+#if DEBUG_COLLAPSE
+  cout<<"Total memory before="<<GetMemUsage()<<endl;
+#endif
+
+  stat.nodes -= collapsedNodes - 1;
+  stat.rayRefs -= totalRayCount - rayCount;
+  
+#if DEBUG_COLLAPSE
+  cout<<"collapsed nodes"<<collapsedNodes<<endl;
+  cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
+  cout<<"Total memory after="<<GetMemUsage()<<endl;
+  cout<<"================================"<<endl;
+#endif
+
+  //  tstat.collapsedNodes += collapsedNodes;
+  //  tstat.collapsedRays += totalRayCount - rayCount;
+    
+  return totalRayCount - rayCount;
+}
+
+
+int
+RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
+{
+  stack<RssTreeNode *> tstack;
+  tstack.push(root);
+
+  Intersectable::NewMail();
+  int pvsSize = 0;
+	
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
+    if (node->IsLeaf()) {
+	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+	  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		  ri != leaf->rays.end();
+		  ri++)
+		if ((*ri).mRay->IsActive()) {
+		  Intersectable *object;
+#if BIDIRECTIONAL_RAY
+		  object = (*ri).mRay->mOriginObject;
+		  if (object && !object->Mailed()) {
+			pvsSize++;
+			object->Mail();
+		  }
+#endif
+		  object = (*ri).mRay->mTerminationObject;
+		  if (object && !object->Mailed()) {
+			pvsSize++;
+			object->Mail();
+		  }
+		}
+	} else {
+	  RssTreeInterior *in = (RssTreeInterior *)node;
+	  if (in->axis < 3) {
+		if (box.Max(in->axis) >= in->position )
+		  tstack.push(in->front);
+				
+		if (box.Min(in->axis) <= in->position )
+		  tstack.push(in->back);
+	  } else {
+		// both nodes for directional splits
+		tstack.push(in->front);
+		tstack.push(in->back);
+	  }
+	}
+  }
+  return pvsSize;
+}
+
+int
+RssTree::CollectPvs(const AxisAlignedBox3 &box,
+					ObjectContainer &pvs
+					) const
+{
+  stack<RssTreeNode *> tstack;
+  tstack.push(root);
+  
+  Intersectable::NewMail();
+  
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
+    if (node->IsLeaf()) {
+	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+	  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+		  ri != leaf->rays.end();
+		  ri++)
+		if ((*ri).mRay->IsActive()) {
+		  Intersectable *object;
+		  object = (*ri).mRay->mTerminationObject;
+		  if (object && !object->Mailed()) {
+			pvs.push_back(object);
+			object->Mail();
+		  }
+		}
+	} else {
+	  RssTreeInterior *in = (RssTreeInterior *)node;
+	  if (in->axis < 3) {
+		if (box.Max(in->axis) >= in->position )
+		  tstack.push(in->front);
+		
+		if (box.Min(in->axis) <= in->position )
+		  tstack.push(in->back);
+	  } else {
+		// both nodes for directional splits
+		tstack.push(in->front);
+		tstack.push(in->back);
+	  }
+	}
+  }
+  return pvs.size();
+}
+
+void
+RssTree::GetTreeStatistics(
+						   float &avgPvsSize,
+						   float &avgRayContribution,
+						   float &avgPvsEntropy,
+						   float &avgRayLengthEntropy,
+						   float &avgImportance
+						   )
+{
+  stack<RssTreeNode *> tstack;
+  tstack.push(root);
+	
+  float sumPvsSize = 0.0f;
+  float sumRayContribution = 0.0f;
+  float sumImportance = 0.0f;
+  float sumPvsEntropy = 0.0f;
+  float sumRayLengthEntropy = 0.0f;
+  
+  int leaves = 0;
+
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+	  leaves++;
+	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+
+	  sumPvsSize += leaf->GetPvsSize();
+	  sumRayContribution += leaf->GetAvgRayContribution();
+	  sumPvsEntropy += leaf->mPvsEntropy;
+	  sumRayLengthEntropy += leaf->mRayLengthEntropy;
+
+	  float imp = leaf->GetImportance();
+	  
+	  if (imp > 1.0f)
+		cout<<"warning imp > 1.0f:"<<imp<<endl;
+	  
+	  sumImportance += imp;
+	  
+
+	} else {
+	  RssTreeInterior *in = (RssTreeInterior *)node;
+	  // both nodes for directional splits
+	  tstack.push(in->front);
+	  tstack.push(in->back);
+	}
+  }
+  
+  avgPvsSize = sumPvsSize/(float)leaves;
+  avgRayContribution = sumRayContribution/(float)leaves;
+  avgPvsEntropy = sumPvsEntropy/(float)leaves;
+  avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
+  avgImportance = sumImportance/(float)leaves;
+}
+
+
+int
+RssTree::GenerateRays(const float ratioPerLeaf,
+					  SimpleRayContainer &rays)
+{
+  stack<RssTreeNode *> tstack;
+  tstack.push(root);
+	
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+	  float c = leaf->GetImportance();
+	  int num = (c*ratioPerLeaf + 0.5);
+	  //			cout<<num<<" ";
+
+	  for (int i=0; i < num; i++) {
+		Vector3 origin = GetBBox(leaf).GetRandomPoint();
+		Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
+		Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
+		//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
+		rays.push_back(SimpleRay(origin, direction));
+	  }
+			
+	} else {
+	  RssTreeInterior *in = (RssTreeInterior *)node;
+	  // both nodes for directional splits
+	  tstack.push(in->front);
+	  tstack.push(in->back);
+	}
+  }
+
+  return rays.size();
+}
+
+void
+RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
+{
+  stack<RssTreeNode *> tstack;
+  tstack.push(root);
+	
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+	  leaves.push_back(leaf);
+	} else {
+	  RssTreeInterior *in = (RssTreeInterior *)node;
+	  // both nodes for directional splits
+	  tstack.push(in->front);
+	  tstack.push(in->back);
+	}
+  }
+}
+
+bool
+RssTree::ValidLeaf(RssTreeLeaf *leaf) const
+{
+  return leaf->rays.size() > termMinRays/4;
+}
+
+
+void
+RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
+						  const int numberOfRays,
+						  SimpleRayContainer &rays)
+{
+  int nrays = leaf->rays.size();
+  for (int i=0; i < numberOfRays; i++) {
+	// pickup 3 random rays
+	int r1 = Random(nrays-1);
+	int r2 = Random(nrays-1);
+	int r3 = Random(nrays-1);
+		
+	Vector3 o1 = leaf->rays[r1].GetOrigin();
+
+	Vector3 o2 = leaf->rays[r2].GetOrigin();
+		
+	Vector3 o3 = leaf->rays[r3].GetOrigin();
+
+	const float overlap = 0.1;
+		
+	Vector3 origin, direction;
+	// generate half of convex combination and half of random rays
+	bool useExtendedConvexCombination = (i > numberOfRays/2);
+	if (useExtendedConvexCombination) {
+	  float w1, w2, w3;
+	  GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
+	  origin = w1*o1 + w2*o2 + w3*o3;
+	  direction =
+		w1*leaf->rays[r1].mRay->GetDir() +
+		w2*leaf->rays[r2].mRay->GetDir() +
+		w3*leaf->rays[r3].mRay->GetDir();
+	  // shift the origin a little bit
+	  origin += direction*0.1f;
+	} else {
+	  origin = GetBBox(leaf).GetRandomPoint();
+	  Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
+	  direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
+	}
+	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
+	rays.push_back(SimpleRay(origin, direction));
+  }
+}
+
+int
+RssTree::GenerateRays(const int numberOfRays,
+					  const int numberOfLeaves,
+					  SimpleRayContainer &rays)
+{
+
+  vector<RssTreeLeaf *> leaves;
+	
+  CollectLeaves(leaves);
+	
+  sort(leaves.begin(),
+	   leaves.end(),
+	   GreaterContribution);
+			 
+
+  float sumContrib = 0.0;
+  int i;
+  int k = 0;
+  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
+	if (ValidLeaf(leaves[i])) {
+	  float c = leaves[i]->GetImportance();
+	  sumContrib += c;
+	  //		cout<<"ray contrib "<<i<<" : "<<c<<endl;
+	  k++;
+	}
+				
+  float avgContrib = sumContrib/numberOfLeaves;
+  float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
+  k = 0;
+  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
+	if (ValidLeaf(leaves[i])) {
+	  k++;
+	  RssTreeLeaf *leaf = leaves[i];
+	  float c = leaf->GetImportance();
+	  int num = (c*ratioPerLeaf + 0.5);
+	  GenerateLeafRays(leaf, num, rays);
+	}
+
+  return rays.size();
+}
+
+
+float
+RssTree::GetAvgPvsSize()
+{
+  stack<RssTreeNode *> tstack;
+  tstack.push(root);
+
+  int sumPvs = 0;
+  int leaves = 0;
+  while (!tstack.empty()) {
+    RssTreeNode *node = tstack.top();
+    tstack.pop();
+		
+    if (node->IsLeaf()) {
+	  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+	  // update pvs size
+	  leaf->UpdatePvsSize();
+	  sumPvs += leaf->GetPvsSize();
+	  leaves++;
+	} else {
+	  RssTreeInterior *in = (RssTreeInterior *)node;
+	  // both nodes for directional splits
+	  tstack.push(in->front);
+	  tstack.push(in->back);
+	}
+  }
+
+
+  return sumPvs/(float)leaves;
+}
+
+	
+float
+RssTreeLeaf::GetImportance()  const
+{
+
+  if (1) {
+	return GetAvgRayContribution();
+	//	return GetPvsSize();
+  } else {
+	// return GetAvgRayContribution()*mEntropyImportance;
+	//return GetAvgRayContribution();
+	return mEntropyImportance;
+  }
+}
+
+
+float
+RssTreeLeaf::ComputePvsEntropy()
+{
+  int samples = 0;
+  Intersectable::NewMail();
+  // set all object as belonging to the fron pvs
+  for(RssTreeNode::RayInfoContainer::iterator ri = rays.begin();
+	  ri != rays.end();
+	  ri++)
+	if ((*ri).mRay->IsActive()) {
+	  Intersectable *object = (*ri).mRay->mTerminationObject;
+	  if (object) {
+		if (!object->Mailed()) {
+		  object->Mail();
+		  object->mCounter = 1;
+		} else
+		  object->mCounter++;
+		samples++;
+	  }
+	}
+  
+  float entropy = 0.0f;
+  
+  if (samples > 1) {
+	Intersectable::NewMail();
+	for(RayInfoContainer::const_iterator ri = rays.begin();
+		ri != rays.end();
+		ri++) 
+	  if ((*ri).mRay->IsActive()) {
+		Intersectable *object = (*ri).mRay->mTerminationObject;
+		if (object) {
+		  if (!object->Mailed()) {
+			object->Mail();
+			float p = object->mCounter/(float)samples;
+			entropy -= p*log(p);
+		  }
+		}
+	  }
+	entropy = entropy/log((float)samples);
+  }
+  else
+	entropy = 1.0f;
+  
+  return entropy;
+}
+
+float
+RssTreeLeaf::ComputeRayLengthEntropy()
+{
+  // get sum of all ray lengths
+  // consider only passing rays or originating rays
+  float sum = 0.0f;
+  int samples = 0;
+  int i=0;
+  for(RayInfoContainer::const_iterator ri = rays.begin();
+      ri != rays.end();
+      ri++) 
+	if ((*ri).mRay->IsActive()) {
+//  	  float s;
+//  	  if (i == 0)
+//  		s = 200;
+//  	  else
+//  		s = 100;
+//  	  i++;
+	  
+	  sum += (*ri).mRay->GetSize();
+	  samples++;
+	}
+  
+  
+  float entropy = 0.0f;
+  
+  if (samples > 1) {
+	i = 0;
+	for(RayInfoContainer::const_iterator ri = rays.begin();
+		ri != rays.end();
+		ri++)
+	  if ((*ri).mRay->IsActive()) {
+//  		float s;
+//  		if (i==0)
+//  		  s = 200;
+//  		else
+//  		  s = 100;
+//  		i++;
+//  		float p = s/sum;
+		float p = (*ri).mRay->GetSize()/sum;
+		entropy -= p*log(p);
+	  }
+	entropy = entropy/log((float)samples);
+  } else
+	entropy = 1.0f;
+  
+  return entropy;
+}
+  
+
+
+void
+RssTreeLeaf::ComputeEntropyImportance()
+{
+  mPvsEntropy = ComputePvsEntropy();
+  mRayLengthEntropy = ComputeRayLengthEntropy();
+  
+  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
+  mEntropyImportance = 1.0f - ComputePvsEntropy();
+  
+  //  cout<<"ei="<<mEntropyImportance<<" ";
+}
+
+
+AxisAlignedBox3
+RssTree::GetShrankedBBox(const RssTreeNode *node) const
+{
+  if (node->parent == NULL)
+	return bbox;
+  
+  if (!node->IsLeaf())
+	return ((RssTreeInterior *)node)->bbox;
+
+  // evaluate bounding box from the ray origins
+  AxisAlignedBox3 box;
+  box.Initialize();
+  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+	  ri != leaf->rays.end();
+	  ri++)
+	if ((*ri).mRay->IsActive()) {
+	  box.Include((*ri).GetOrigin());
+	}
+  return box;
+}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp	(revision 447)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp	(revision 447)
@@ -0,0 +1,665 @@
+#include "SceneGraph.h"
+#include "KdTree.h"
+#include "RssPreprocessor.h"
+#include "X3dExporter.h"
+#include "Environment.h"
+#include "MutualVisibility.h"
+#include "Polygon3.h"
+#include "ViewCell.h"
+#include "VssRay.h"
+#include "RssTree.h"
+
+
+static bool useViewSpaceBox = true;//true;
+static bool use2dSampling = false;
+static bool fromBoxVisibility = true;
+
+static bool exportPvs = false;
+
+RssPreprocessor::RssPreprocessor():
+  mPass(0),
+  mVssRays()
+{
+  // this should increase coherence of the samples
+  environment->GetIntValue("RssPreprocessor.samplesPerPass", mSamplesPerPass);
+  environment->GetIntValue("RssPreprocessor.initialSamples", mInitialSamples);
+  environment->GetIntValue("RssPreprocessor.vssSamples", mRssSamples);
+  environment->GetIntValue("RssPreprocessor.vssSamplesPerPass", mRssSamplesPerPass);
+  environment->GetBoolValue("RssPreprocessor.useImportanceSampling", mUseImportanceSampling);
+  environment->GetIntValue("BspTree.Construction.samples", mBspConstructionSamples);
+	
+  mStats.open("stats.log");
+}
+
+RssPreprocessor::~RssPreprocessor()
+{
+  CLEAR_CONTAINER(mVssRays);
+}
+
+void
+RssPreprocessor::SetupRay(Ray &ray, 
+						  const Vector3 &point, 
+						  const Vector3 &direction
+						  )
+{
+  ray.intersections.clear();
+  // do not store anything else then intersections at the ray
+  ray.Init(point, direction, Ray::LOCAL_RAY);
+}
+
+int
+RssPreprocessor::CastRay(
+						 Vector3 &viewPoint,
+						 Vector3 &direction,
+						 VssRayContainer &vssRays
+						 )
+{
+  int hits = 0;
+  static Ray ray;
+  AxisAlignedBox3 box = mKdTree->GetBox();
+
+  AxisAlignedBox3 sbox = box;
+  sbox.Enlarge(Vector3(-Limits::Small));
+  if (!sbox.IsInside(viewPoint))
+	return 0;
+	
+  SetupRay(ray, viewPoint, direction);
+  // cast ray to KD tree to find intersection with other objects
+  Intersectable *objectA, *objectB;
+  Vector3 pointA, pointB;
+  float bsize = Magnitude(box.Size());
+  
+  if (mKdTree->CastRay(ray)) {
+	objectA = ray.intersections[0].mObject;
+	pointA = ray.Extrap(ray.intersections[0].mT);
+  } else {
+	objectA = NULL;
+	// compute intersection with the scene bounding box
+	float tmin, tmax;
+	box.ComputeMinMaxT(ray, &tmin, &tmax);
+	if (tmax > bsize) {
+	  //cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
+	  //cerr<<"ray"<<ray<<endl;
+	}
+	pointA = ray.Extrap(tmax);
+  }
+
+  bool detectEmptyViewSpace = true;
+	
+  if (detectEmptyViewSpace) {
+	SetupRay(ray, pointA, -direction);
+  } else
+	SetupRay(ray, viewPoint, -direction);
+	
+	
+  if (mKdTree->CastRay(ray)) {
+	objectB = ray.intersections[0].mObject;
+	pointB = ray.Extrap(ray.intersections[0].mT);
+  } else {
+	objectB = NULL;
+	float tmin, tmax;
+	box.ComputeMinMaxT(ray, &tmin, &tmax);
+	if (tmax > bsize) {
+	  //cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
+	  //cerr<<"ray"<<ray<<endl;
+	}
+	pointB = ray.Extrap(tmax);
+  }
+
+  VssRay *vssRay  = NULL;
+
+  bool validSample = true;
+  if (detectEmptyViewSpace) {
+	// check if the viewpoint lies on the line segment AB
+	if (Distance(pointA, pointB) <
+		Distance(viewPoint, pointA) + Distance(viewPoint, pointB) - Limits::Small) {
+	  validSample = false;
+	}
+  }
+	
+  if (validSample) {
+	if (objectA) {
+	  vssRay = new VssRay(pointB,
+						  pointA,
+						  objectB,
+						  objectA);
+	  vssRays.push_back(vssRay);
+	  hits ++;
+	}
+	
+	if (objectB) {
+	  vssRay = new VssRay(pointA,
+						  pointB,
+						  objectA,
+						  objectB);
+	  vssRays.push_back(vssRay);
+	  hits ++;
+	}
+  }
+	
+  return hits;
+}
+
+
+Vector3
+RssPreprocessor::GetViewpoint(AxisAlignedBox3 *viewSpaceBox)
+{
+  AxisAlignedBox3 box;
+	
+  if (viewSpaceBox)
+	box =*viewSpaceBox;
+  else 
+	box = mKdTree->GetBox();
+	
+  // shrink the box in the y direction
+  return box.GetRandomPoint();
+}
+
+Vector3
+RssPreprocessor::GetDirection(const Vector3 &viewpoint,
+							  AxisAlignedBox3 *viewSpaceBox
+							  )
+{
+  Vector3 point;
+  if (!use2dSampling) {
+	Vector3 normal;
+	int i = RandomValue(0, mObjects.size()-1);
+	Intersectable *object = mObjects[i];
+	object->GetRandomSurfacePoint(point, normal);
+  } else {
+	AxisAlignedBox3 box;
+		
+	if (viewSpaceBox)
+	  box =*viewSpaceBox;
+	else 
+	  box = mKdTree->GetBox();
+		
+	point = box.GetRandomPoint();
+	point.y = viewpoint.y;
+  }
+	
+  return point - viewpoint;
+}
+
+int
+RssPreprocessor::GenerateImportanceRays(RssTree *rssTree,
+										const int desiredSamples,
+										SimpleRayContainer &rays
+										)
+{
+  int num;
+  
+  float avgPvsSize;
+  float avgRayContribution;
+  float avgPvsEntropy;
+  float avgRayLengthEntropy;
+  float avgImportance;
+  
+  rssTree->GetTreeStatistics(
+							 avgPvsSize,
+							 avgRayContribution,
+							 avgPvsEntropy,
+							 avgRayLengthEntropy,
+							 avgImportance);
+  
+  cout<<
+	"#AVG_PVS_SIZE\n"<<avgPvsSize<<endl<<
+	"#AVG_RAY_CONTRIB\n"<<avgRayContribution<<endl<<
+	"#AVG_PVS_ENTROPY\n"<<avgPvsEntropy<<endl<<
+	"#AVG_RAY_LENGTH_ENTROPY\n"<<avgRayLengthEntropy<<endl<<
+	"#AVG_IMPORTANCE\n"<<avgImportance<<endl;
+  
+  if (0) {
+	float p = desiredSamples/(float)(avgRayContribution*rssTree->stat.Leaves());
+	num = rssTree->GenerateRays(p, rays);
+  } else {
+	int leaves = rssTree->stat.Leaves()/2;
+	num = rssTree->GenerateRays(desiredSamples, leaves, rays);
+  }
+	
+  cout<<"Generated "<<num<<" rays."<<endl;
+	
+  return num;
+}
+
+
+bool
+RssPreprocessor::ExportRays(const char *filename,
+							const VssRayContainer &vssRays,
+							const int number
+							)
+{
+  cout<<"Exporting vss rays..."<<endl<<flush;
+	
+  float prob = number/(float)vssRays.size();
+
+
+  Exporter *exporter = NULL;
+  exporter = Exporter::GetExporter(filename);
+  //	exporter->SetWireframe();
+  //	exporter->ExportKdTree(*mKdTree);
+  exporter->SetFilled();
+  exporter->ExportScene(mSceneGraph->mRoot);
+  exporter->SetWireframe();
+
+  if (mViewSpaceBox) {
+	exporter->SetForcedMaterial(RgbColor(1,0,1));
+	exporter->ExportBox(*mViewSpaceBox);
+	exporter->ResetForcedMaterial();
+  }
+	
+  VssRayContainer rays;	for (int i=0; i < vssRays.size(); i++)
+	if (RandomValue(0,1) < prob)
+	  rays.push_back(vssRays[i]);
+
+  exporter->ExportRays(rays, RgbColor(1, 0, 0));
+	
+  delete exporter;
+
+  cout<<"done."<<endl<<flush;
+
+  return true;
+}
+
+
+bool
+RssPreprocessor::ExportRssTree(char *filename,
+							   RssTree *tree,
+							   const Vector3 &dir
+							   )
+{
+  Exporter *exporter = Exporter::GetExporter(filename);
+  exporter->SetFilled();
+  exporter->ExportScene(mSceneGraph->mRoot);
+  //  exporter->SetWireframe();
+  bool result = exporter->ExportRssTree2( *tree, dir );
+  delete exporter;
+  return result;
+}
+
+bool
+RssPreprocessor::ExportRssTreeLeaf(char *filename,
+								   RssTree *tree,
+								   RssTreeLeaf *leaf)
+{
+  Exporter *exporter = NULL;
+  exporter = Exporter::GetExporter(filename);
+  exporter->SetWireframe();
+  exporter->ExportKdTree(*mKdTree);
+	
+  if (mViewSpaceBox) {
+	exporter->SetForcedMaterial(RgbColor(1,0,0));
+	exporter->ExportBox(*mViewSpaceBox);
+	exporter->ResetForcedMaterial();
+  }
+	
+  exporter->SetForcedMaterial(RgbColor(0,0,1));
+  exporter->ExportBox(tree->GetBBox(leaf));
+  exporter->ResetForcedMaterial();
+	
+  VssRayContainer rays[4];
+  for (int i=0; i < leaf->rays.size(); i++) {
+	int k = leaf->rays[i].GetRayClass();
+	rays[k].push_back(leaf->rays[i].mRay);
+  }
+	
+  // SOURCE RAY
+  exporter->ExportRays(rays[0], RgbColor(1, 0, 0));
+  // TERMINATION RAY
+  exporter->ExportRays(rays[1], RgbColor(1, 1, 1));
+  // PASSING_RAY
+  exporter->ExportRays(rays[2], RgbColor(1, 1, 0));
+  // CONTAINED_RAY
+  exporter->ExportRays(rays[3], RgbColor(0, 0, 1));
+
+  delete exporter;
+  return true;
+}
+
+void
+RssPreprocessor::ExportRssTreeLeaves(RssTree *tree, const int number)
+{
+  vector<RssTreeLeaf *> leaves;
+  tree->CollectLeaves(leaves);
+
+  int num = 0;
+  int i;
+  float p = number / (float)leaves.size();
+  for (i=0; i < leaves.size(); i++) {
+	if (RandomValue(0,1) < p) {
+	  char filename[64];
+	  sprintf(filename, "rss-leaf-%04d.x3d", num);
+	  ExportRssTreeLeaf(filename, tree, leaves[i]);
+	  num++;
+	}
+	if (num >= number)
+	  break;
+  }
+}
+
+
+float
+RssPreprocessor::GetAvgPvsSize(RssTree *tree,
+							   const vector<AxisAlignedBox3> &viewcells
+							   )
+{
+  vector<AxisAlignedBox3>::const_iterator it, it_end = viewcells.end();
+
+  int sum = 0;
+  for (it = viewcells.begin(); it != it_end; ++ it)
+	sum += tree->GetPvsSize(*it);
+	
+  return sum/(float)viewcells.size();
+}
+
+
+void
+RssPreprocessor::ExportPvs(char *filename,
+						   RssTree *rssTree
+						   )
+{
+  ObjectContainer pvs;
+  if (rssTree->CollectRootPvs(pvs)) {
+	Exporter *exporter = Exporter::GetExporter(filename);
+	exporter->SetFilled();
+	exporter->ExportGeometry(pvs);
+	delete exporter;
+  }
+}
+
+bool
+RssPreprocessor::ComputeVisibility()
+{
+  
+  mSceneGraph->CollectObjects(&mObjects);
+	
+  long startTime = GetTime();
+  
+  int totalSamples = 0;
+
+
+  AxisAlignedBox3 *box = new AxisAlignedBox3(mKdTree->GetBox());
+
+  if (fromBoxVisibility) {
+	float size = 0.05f;
+	float s = 0.5f - size;
+	float olds = Magnitude(box->Size());
+	box->Enlarge(box->Size()*Vector3(-s));
+	Vector3 translation = Vector3(-olds*0.1f, 0, 0);
+	box->SetMin(box->Min() + translation);
+	box->SetMax(box->Max() + translation);
+  } else {
+		
+	// sample city like heights
+	box->SetMin(1, box->Min(1) + box->Size(1)*0.2);
+	box->SetMax(1, box->Min(1) + box->Size(1)*0.3);
+  }
+
+  if (use2dSampling)
+	box->SetMax(1, box->Min(1));
+	
+  if (useViewSpaceBox)
+	mViewSpaceBox = box;
+  else
+	mViewSpaceBox = NULL;
+		
+
+  RssTree *rssTree = NULL;
+
+  RayContainer bspRays;
+
+  while (totalSamples < mInitialSamples) {
+	int passContributingSamples = 0;
+	int passSampleContributions = 0;
+	int passSamples = 0;
+	int index = 0;
+		
+	int sampleContributions;
+		
+	int s = Min(mSamplesPerPass, mInitialSamples);
+	for (int k=0; k < s; k++) {
+			
+	  Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
+	  Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
+			
+	  sampleContributions = CastRay(viewpoint, direction, mVssRays);
+			
+			
+	  //-- CORR matt: put block inside loop
+	  if (sampleContributions) {
+		passContributingSamples ++;
+		passSampleContributions += sampleContributions;
+	  }
+	  passSamples++;
+	  totalSamples++;
+	}
+    
+	mPass++;
+		
+	int pvsSize = 0;
+	float avgRayContrib = (passContributingSamples > 0) ? 
+	  passSampleContributions/(float)passContributingSamples : 0;
+		
+	cout << "#Pass " << mPass << " : t = " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+	cout << "#TotalSamples=" << totalSamples/1000 
+		 << "k   #SampleContributions=" << passSampleContributions << " (" 
+		 << 100*passContributingSamples/(float)passSamples<<"%)" << " avgPVS="
+		 << pvsSize/(float)mObjects.size() << endl 
+		 << "avg ray contrib=" << avgRayContrib << endl;
+		
+	mStats <<
+	  "#Pass\n" <<mPass<<endl<<
+	  "#Time\n" << TimeDiff(startTime, GetTime())*1e-3 << endl<<
+	  "#TotalSamples\n" << totalSamples<< endl<<
+	  "#SampleContributions\n" << passSampleContributions << endl << 
+	  "#PContributingSamples\n"<<100*passContributingSamples/(float)passSamples<<endl <<
+	  "#AvgPVS\n"<< pvsSize/(float)mObjects.size() << endl <<
+	  "#AvgRayContrib\n" << avgRayContrib << endl;
+  }
+  
+  cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
+  cout << "#totalRayStackSize=" << mVssRays.size() << endl <<flush;
+  
+  int numExportRays = 10000;
+  //  int numExportRays = 0;
+
+  if (numExportRays) {
+	char filename[64];
+	sprintf(filename, "rss-rays-initial.x3d");
+	ExportRays(filename, mVssRays, numExportRays);
+  }
+  
+  if (exportPvs) {
+	ExportPvs("rss-pvs-initial.x3d", rssTree);
+  }
+  
+  //-- 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);
+	
+  if (fromBoxVisibility)
+	rssTree->Construct(mVssRays, mViewSpaceBox);
+  else
+	rssTree->Construct(mVssRays, NULL);
+
+  cout<<"RssTree root PVS size = "<<rssTree->GetRootPvsSize()<<endl;
+  
+  ExportRssTree("rss-tree-100.x3d", rssTree, Vector3(1,0,0));
+  ExportRssTree("rss-tree-001.x3d", rssTree, Vector3(0,0,1));
+  ExportRssTree("rss-tree-101.x3d", rssTree, Vector3(1,0,1));
+  ExportRssTree("rss-tree-101m.x3d", rssTree, Vector3(-1,0,-1));
+
+  ExportRssTreeLeaves(rssTree, 10);
+
+  // viewcells->UpdatePVS(newVssRays);
+  // get viewcells as kd tree boxes
+  vector<AxisAlignedBox3> kdViewcells;
+  if (0) {
+	vector<KdLeaf *> leaves;
+	mKdTree->CollectLeaves(leaves);
+	vector<KdLeaf *>::const_iterator it;
+	int targetLeaves = 50;
+	float prob = targetLeaves/(float)leaves.size();
+	for (it = leaves.begin(); it != leaves.end(); ++it)
+	  if (RandomValue(0.0f,1.0f) < prob)
+		kdViewcells.push_back(mKdTree->GetBox(*it));
+		
+	float avgPvs = GetAvgPvsSize(rssTree, kdViewcells);
+	cout<<"Initial average PVS size = "<<avgPvs<<endl;
+  }
+
+	
+  int samples = 0;
+  int pass = 0;
+  while (1) {
+	int num = mRssSamplesPerPass;
+	SimpleRayContainer rays;
+	VssRayContainer vssRays;
+		
+	if (!mUseImportanceSampling) {
+	  for (int j=0; j < num; j++) {
+		Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
+		Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
+		rays.push_back(SimpleRay(viewpoint, direction));
+	  }
+	} else {
+	  num = GenerateImportanceRays(rssTree, num, rays);
+	}
+		
+		
+	for (int i=0; i < rays.size(); i++)
+	  CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays);
+		
+	rssTree->AddRays(vssRays);
+		
+	if (1) {
+	  int subdivided = rssTree->UpdateSubdivision();
+	  cout<<"subdivided leafs = "<<subdivided<<endl;
+	}
+
+	float avgPvs = GetAvgPvsSize(rssTree, kdViewcells);
+	cout<<"Average PVS size = "<<avgPvs<<endl;
+
+	if (numExportRays) {
+	  char filename[64];
+	  if (mUseImportanceSampling)
+		sprintf(filename, "rss-rays-i%04d.x3d", pass);
+	  else
+		sprintf(filename, "rss-rays-%04d.x3d", pass);
+	  
+	  ExportRays(filename, vssRays, numExportRays);
+	}
+
+	// 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();
+	cout<<"*****************************\n";
+	cout<<samples<<" avgPVS ="<<pvs<<endl;
+	cout<<"RssTree root PVS size = "<<rssTree->GetRootPvsSize()<<endl;
+	cout<<"*****************************\n";
+
+	if (exportPvs) {
+	  char filename[64];
+	  sprintf(filename, "rss-pvs-%04d.x3d", pass);
+	  ExportPvs(filename, rssTree);
+	}
+
+	if (samples >= mRssSamples)
+	  break;
+
+	
+	pass++;
+  }
+  
+  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;
+}
+
+void
+RssPreprocessor::CastRay(const BspTree &tree, const VssRay & vssRay)
+{
+  //-- cast ray to BSP tree to get intersection with view cells
+  Ray ray(vssRay);
+  mBspTree->CastRay(ray);
+				
+  //if (0 && ray.sourceObject.mObject)
+  //sampleContributions += 
+  //AddObjectSamples(ray.sourceObject.mObject, ray);
+			
+  if (!ray.intersections.empty()) // second intersection found
+  {
+	  //sampleContributions += 
+	  AddObjectSamples(ray.intersections[0].mObject, ray);
+  }
+}
+
+int
+RssPreprocessor::AddObjectSamples(Intersectable *obj, const Ray &ray)
+{
+  int contributingSamples = 0;
+  int j;
+ 
+  // object can be seen from the view cell => add to view cell pvs
+  for (j=0; j < ray.bspIntersections.size(); ++ j)
+	{ 	
+	  BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
+	  // if ray not in unbounded space
+	  if (leaf->GetViewCell() != &mUnbounded)
+		contributingSamples += 
+		  leaf->GetViewCell()->GetPvs().AddSample(obj);
+	}
+ 
+  // rays passing through this viewcell
+  if (mPass > 1)
+	for (j=1; j < ((int)ray.bspIntersections.size() - 1); ++ j) 
+	  {
+		BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
+		
+		if (leaf->GetViewCell() != &mUnbounded)
+		  leaf->GetViewCell()->
+			AddPassingRay(ray, contributingSamples ? 1 : 0);
+	  }
+  
+  return contributingSamples;
+}
Index: /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.h
===================================================================
--- /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.h	(revision 447)
+++ /trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.h	(revision 447)
@@ -0,0 +1,113 @@
+#ifndef _RssPreprocessor_H__
+#define _RssPreprocessor_H__
+
+#include <fstream>
+using namespace std;
+
+#include "Preprocessor.h"
+#include "VssRay.h"
+class RssTree;
+class RssTreeLeaf;
+
+/** Sampling based visibility preprocessing. The implementation is based on heuristical
+    sampling of view space */
+class RssPreprocessor : public Preprocessor {
+public:
+  int mPass;
+  int mSamplesPerPass;
+  int mRssSamplesPerPass;
+  int mInitialSamples;
+  int mRssSamples;
+  bool mUseImportanceSampling;
+
+  AxisAlignedBox3 *mViewSpaceBox;
+
+  ofstream mStats;
+	
+  ObjectContainer mObjects;
+
+  // rays cast during the processing
+  VssRayContainer mVssRays;
+	
+  RssPreprocessor();
+  ~RssPreprocessor();
+
+  virtual bool ComputeVisibility();
+
+  Vector3
+  GetViewpoint(AxisAlignedBox3 *viewSpaceBox);
+
+  Vector3
+  GetDirection(const Vector3 &viewpoint,
+			   AxisAlignedBox3 *viewSpaceBox
+			   );
+	
+  void
+  SetupRay(Ray &ray, 
+		   const Vector3 &point, 
+		   const Vector3 &direction
+		   );
+
+
+	
+  int
+  CastRay(
+		  Vector3 &viewPoint,
+		  Vector3 &direction,
+		  VssRayContainer &vssRays
+
+		  );
+		
+
+	
+  virtual bool BuildBspTree() { return false; }
+
+
+  bool
+  ExportRays(const char *filename,
+			 const VssRayContainer &vssRays,
+			 const int number
+			 );
+
+  int
+  GenerateImportanceRays(RssTree *vssTree,
+						 const int desiredSamples,
+						 SimpleRayContainer &rays
+						 );
+
+
+  bool
+  ExportRssTreeLeaf(char *filename,
+					RssTree *tree,
+					RssTreeLeaf *leaf);
+
+  void
+  ExportRssTreeLeaves(RssTree *tree, const int number);
+	
+
+  bool
+  ExportRssTree(char *filename,
+				RssTree *tree,
+				const Vector3 &dir
+				);
+
+
+  float
+  GetAvgPvsSize(RssTree *tree,
+				const vector<AxisAlignedBox3> &viewcells
+				);
+  
+  void CastRay(const BspTree &tree, const VssRay & vssRay);
+	
+  int AddObjectSamples(Intersectable *obj, const Ray &ray);
+
+  void
+  ExportPvs(char *filename,
+			RssTree *rssTree
+			);
+	
+};
+
+
+
+#endif
