Changeset 408 for trunk/VUT

Timestamp:

11/14/05 11:12:38 (19 years ago)

Author:

mattausch

Message:

 

Location:

trunk/VUT/GtpVisibilityPreprocessor

Files:

• scripts/Preprocessor.vcproj (modified) (1 diff)
• src/Environment.cpp (modified) (2 diffs)
• src/VspKdTree.cpp (moved) (moved from trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.cpp) (2 diffs)
• src/VspKdTree.h (moved) (moved from trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellKd.h) (5 diffs)

trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj

@@ -291 +291 @@
                         </File>
                         <File
-                                RelativePath="..\src\ViewCellKd.cpp">
-                        </File>
-                        <File
-                                RelativePath="..\src\ViewCellKd.h">
+                                RelativePath="..\src\VspKdTree.cpp">
+                        </File>
+                        <File
+                                RelativePath="..\src\VspKdTree.h">
                         </File>
                         <File

trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp

@@ -1322 +1322 @@
 
 
-        RegisterOption("VssTree.maxDepth", optInt, "kd_depth=", "12");
+        RegisterOption("VspKdTree.maxDepth", optInt, "kd_depth=", "12");
+  RegisterOption("VspKdTree.minPvs", optInt, "kd_minpvs=", "1");
+  RegisterOption("VspKdTree.minRays", optInt, "kd_minrays=", "10");
+        RegisterOption("VspKdTree.maxCostRatio", optFloat, "maxcost=", "0.95");
+        RegisterOption("VspKdTree.maxRayContribution", optFloat, "maxraycontrib=", "0.5");
+
+        RegisterOption("VspKdTree.epsilon", optFloat, "kd_eps=", "1e-6");
+  RegisterOption("VspKdTree.ct_div_ci", optFloat, "kd_ctdivci=", "1.0");
+  RegisterOption("VspKdTree.randomize", optBool, "randomize", "false");
+  RegisterOption("VspKdTree.splitType", optString, "split=", "queries");
+  RegisterOption("VspKdTree.numberOfEndPointDomains", optInt, "endpoints=", "10000");
+
+  RegisterOption("VspKdTree.minSize", optFloat, "minsize=", "0.001");
+
+  RegisterOption("VspKdTree.maxTotalMemory", optFloat, "mem=", "60.0");
+  RegisterOption("VspKdTree.maxStaticMemory", optFloat, "statmem=", "8.0");
+
+  RegisterOption("VspKdTree.queryType", optString, "qtype=", "static");
+
+  RegisterOption("VspKdTree.queryPosWeight", optFloat, "qposweight=", "0.0");
+  RegisterOption("VspKdTree.useRefDirSplits", optBool, "refdir", "false");
+  RegisterOption("VspKdTree.refDirAngle", optFloat, "refangle=", "10");
+  RegisterOption("VspKdTree.refDirBoxMaxSize", optFloat, "refboxsize=", "0.1");
+  RegisterOption("VspKdTree.accessTimeThreshold", optInt, "accesstime=", "1000");
+  RegisterOption("VspKdTree.minCollapseDepth", optInt, "colldepth=", "4");
+
+
+
+
+
+
+
+
+
+
+        RegisterOption("VssTree.maxDepth", optInt, "kd_depth=", "12");
   RegisterOption("VssTree.minPvs", optInt, "kd_minpvs=", "1");
   RegisterOption("VssTree.minRays", optInt, "kd_minrays=", "10");
@@ -1347 +1382 @@
   RegisterOption("VssTree.accessTimeThreshold", optInt, "accesstime=", "1000");
   RegisterOption("VssTree.minCollapseDepth", optInt, "colldepth=", "4");
-
 }
 

trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp

@@ +1 @@
+
+// ================================================================
+// $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 <queue>
+#include <fstream>
+#include <string>
+
+#include "VspKdTree.h"
+
 #include "Environment.h"
-#include "Mesh.h"
-#include "ViewCellKd.h"
-
-
-int ViewCellKdNode::mailID = 1;
-
-
-ViewCellKdNode::ViewCellKdNode(ViewCellKdInterior *parent):mParent(parent), mailbox(0)
-{
-  if (parent)
-    mDepth = parent->mDepth+1;
+#include "VssRay.h"
+#include "Intersectable.h"
+#include "Ray.h"
+
+
+#define DEBUG_DIR_SPLIT 0
+
+
+
+// Static variables
+int
+VspKdTreeLeaf::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
+VspKdTree::VspKdTree()
+{
+  environment->GetIntValue("VspKdTree.maxDepth", termMaxDepth);
+        environment->GetIntValue("VspKdTree.minPvs", termMinPvs);
+        environment->GetIntValue("VspKdTree.minRays", termMinRays);
+        environment->GetFloatValue("VspKdTree.maxRayContribution", termMaxRayContribution);
+  environment->GetFloatValue("VspKdTree.maxCostRatio", termMaxCostRatio);
+
+  environment->GetFloatValue("VspKdTree.minSize", termMinSize);
+  termMinSize = sqr(termMinSize);
+        
+  environment->GetFloatValue("VspKdTree.refDirBoxMaxSize", refDirBoxMaxSize);
+  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
+  
+  environment->GetFloatValue("VspKdTree.epsilon", epsilon);
+  environment->GetFloatValue("VspKdTree.ct_div_ci", ct_div_ci);
+        
+  environment->GetFloatValue("VspKdTree.maxTotalMemory", maxTotalMemory);
+  environment->GetFloatValue("VspKdTree.maxStaticMemory", maxStaticMemory);
+  
+  
+
+
+  float refDirAngle;
+  environment->GetFloatValue("VspKdTree.refDirAngle", refDirAngle);
+
+  environment->GetIntValue("VspKdTree.accessTimeThreshold", accessTimeThreshold);
+  //= 1000;
+  environment->GetIntValue("VspKdTree.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("VspKdTree.splitType", sname);
+  string name(sname);
+        
+  if (name.compare("regular") == 0)
+    splitType = ESplitRegular;
   else
-    mDepth = 0;
-}
-
-ViewCellKd::ViewCellKd()
-{
-  mRoot = new ViewCellKdLeaf(NULL, 0);
-  environment->GetIntValue("ViewCellKd.Termination.maxDepth", mTermMaxDepth);
-  environment->GetIntValue("ViewCellKd.Termination.minCost", mTermMinCost);
-  environment->GetFloatValue("ViewCellKd.Termination.maxCostRatio", mMaxCostRatio);
-  environment->GetFloatValue("ViewCellKd.Termination.ct_div_ci", mCt_div_ci);
-  environment->GetFloatValue("ViewCellKd.splitBorder", mSplitBorder);
-
-  environment->GetBoolValue("ViewCellKd.sahUseFaces", mSahUseFaces);
-
-  char splitType[64];
-  environment->GetStringValue("ViewCellKd.splitMethod", splitType);
-  
-  mSplitMethod = SPLIT_SPATIAL_MEDIAN;
-  if (strcmp(splitType, "spatialMedian") == 0)
-    mSplitMethod = SPLIT_SPATIAL_MEDIAN;
-  else
-    if (strcmp(splitType, "objectMedian") == 0)
-      mSplitMethod = SPLIT_OBJECT_MEDIAN;
-    else
-      if (strcmp(splitType, "SAH") == 0)
-        mSplitMethod = SPLIT_SAH;
-      else {
-        cerr<<"Wrong kd split type "<<splitType<<endl;
-        exit(1);
-      }
-  splitCandidates = NULL;
-}
-
-bool
-ViewCellKd::Construct()
-{
-  if (!splitCandidates)
+    if (name.compare("heuristic") == 0)
+      splitType = ESplitHeuristic;
+                else {
+                        cerr<<"Invalid VspKdTree split type "<<name<<endl;
+                        exit(1);
+                }
+        
+  environment->GetBoolValue("VspKdTree.randomize", randomize);
+
+  root = NULL;
+  
+  splitCandidates = new vector<SortableEntry>;
+}
+
+
+VspKdTree::~VspKdTree()
+{
+  if (root)
+    delete root;
+}
+
+
+
+
+void
+VspKdStatistics::Print(ostream &app) const
+{
+  app << "===== VspKdTree 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_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 VspKdTree statistics ==========\n";
+
+}
+
+
+void
+VspKdTreeLeaf::UpdatePvsSize()
+{
+        if (!mValidPvs) {
+                Intersectable::NewMail();
+                int pvsSize = 0;
+                for(VspKdTreeNode::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();
+                                }
+                        }
+                mPvsSize = pvsSize;
+                mValidPvs = true;
+        }
+}
+
+
+void
+VspKdTree::Construct(
+                                                                         VssRayContainer &rays,
+                                                                         AxisAlignedBox3 *forcedBoundingBox
+                                                                         )
+{
+  stat.Start();
+  
+        maxMemory = maxStaticMemory;
+
+  if (root)
+    delete root;
+
+  root = new VspKdTreeLeaf(NULL, rays.size());
+  // first construct a leaf that will get subdivide
+  VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) root;
+
+  stat.nodes = 1;
+  
+  bbox.Initialize();
+  dirBBox.Initialize();
+  
+  for(VssRayContainer::const_iterator ri = rays.begin();
+      ri != rays.end();
+      ri++) {
+    leaf->AddRay(VspKdTreeNode::RayInfo(*ri));
+
+    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();
+  stat.initialPvsSize = leaf->GetPvsSize();
+  // Subdivide();
+  root = Subdivide(TraversalData(leaf, bbox, 0));
+
+  if (splitCandidates) {
+    // force realease of this vector
+    delete splitCandidates;
     splitCandidates = new vector<SortableEntry>;
-
-  // first construct a leaf that will get subdivide
-  ViewCellKdLeaf *leaf = (ViewCellKdLeaf *) mRoot;
-
-  mStat.nodes = 1;
-
-  mBox.Initialize();
-  
-  ObjectContainer::const_iterator mi;
-  for ( mi = leaf->mObjects.begin();
-        mi != leaf->mObjects.end();
-        mi++) {
-    mBox.Include((*mi)->GetBox());
-  }
-
-  cout <<"ViewCellKd Root Box:"<< mBox<<endl;
-  mRoot = Subdivide(TraversalData(leaf, mBox, 0));
-
-  // remove the allocated array
-  delete splitCandidates;
-
-  return true;
-}
-
-ViewCellKdNode *
-ViewCellKd::Subdivide(const TraversalData &tdata)
-{
-
-  ViewCellKdNode *result = NULL;
+  }
+  
+  stat.Stop();
+
+        stat.Print(cout);
+        cout<<"#Total memory="<<GetMemUsage()<<endl;
+
+}
+
+
+
+VspKdTreeNode *
+VspKdTree::Subdivide(const TraversalData &tdata)
+{
+  VspKdTreeNode *result = NULL;
 
   priority_queue<TraversalData> tStack;
-  //  stack<STraversalData> tStack;
+  //  stack<TraversalData> tStack;
   
   tStack.push(tdata);
-  AxisAlignedBox3 backBox, frontBox;
-
-  
+
+  AxisAlignedBox3 backBox;
+  AxisAlignedBox3 frontBox;
+
+  
+        int lastMem = 0;
   while (!tStack.empty()) {
 
-#if 0
-    if ( GetMemUsage() > maxMemory ) {
+                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();
+                                EvaluateLeafStats(tStack.top());
+                                tStack.pop();
       }
       break;
     }
-#endif
-
+    
     TraversalData data = tStack.top();
     tStack.pop();
     
-    ViewCellKdNode *node = SubdivideNode((ViewCellKdLeaf *) data.mNode,
-                                 data.mBox,
-                                 backBox,
-                                 frontBox
-                                 );
+    VspKdTreeNode *node = SubdivideNode((VspKdTreeLeaf *) data.node,
+                                                                                                                                                        data.bbox,
+                                                                                                                                                        backBox,
+                                                                                                                                                        frontBox
+                                                                                                                                                        );
     if (result == NULL)
       result = node;
     
     if (!node->IsLeaf()) {
-
-      ViewCellKdInterior *interior = (ViewCellKdInterior *) node;
+                        
+      VspKdTreeInterior *interior = (VspKdTreeInterior *) node;
       // push the children on the stack
-      tStack.push(TraversalData(interior->mBack, backBox, data.mDepth+1));
-      tStack.push(TraversalData(interior->mFront, frontBox, data.mDepth+1));
-      
+      tStack.push(TraversalData(interior->back, backBox, data.depth+1));
+      tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
+                        
     } else {
       EvaluateLeafStats(data);
     }
   }
-  
+
   return result;
-
-}
-
-
-
-bool
-ViewCellKd::TerminationCriteriaMet(const ViewCellKdLeaf *leaf)
-{
-  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
-  return
-    (leaf->mObjects.size() <= mTermMinCost) ||
-    (leaf->mDepth >= mTermMaxDepth);
-  
-}
-
-
+}
+
+
+// returns selected plane for subdivision
 int
-ViewCellKd::SelectPlane(ViewCellKdLeaf *leaf,
-                    const AxisAlignedBox3 &box,
-                    float &position
-                    )
-{
-  int axis = -1;
-  
-  switch (mSplitMethod)
-    {
-    case SPLIT_SPATIAL_MEDIAN: {
-      axis = box.Size().DrivingAxis();
-      position = (box.Min()[axis] + box.Max()[axis])*0.5f;
-      break;
+VspKdTree::SelectPlane(
+                                                                                 VspKdTreeLeaf *leaf,
+                                                                                 const AxisAlignedBox3 &box,
+                                                                                 float &position,
+                                                                                 int &raysBack,
+                                                                                 int &raysFront,
+                                                                                 int &pvsBack,
+                                                                                 int &pvsFront
+                                                                                 )
+{
+
+  int minDirDepth = 6;
+  int axis;
+        float costRatio;
+        
+  if (splitType == ESplitRegular) {
+                costRatio = BestCostRatioRegular(leaf,
+                                                                                                                                                 axis,
+                                                                                                                                                 position,
+                                                                                                                                                 raysBack,
+                                                                                                                                                 raysFront,
+                                                                                                                                                 pvsBack,
+                                                                                                                                                 pvsFront
+                                                                                                                                                 );
+
+        } else {
+    if (splitType == ESplitHeuristic)
+      costRatio = BestCostRatioHeuristic(leaf,
+                                                                                                                                                                 axis,
+                                                                                                                                                                 position,
+                                                                                                                                                                 raysBack,
+                                                                                                                                                                 raysFront,
+                                                                                                                                                                 pvsBack,
+                                                                                                                                                                 pvsFront
+                                                                                                                                                                 );
+                else {
+                        cerr<<"VspKdTree: Unknown split heuristics\n";
+                        exit(1);
+                }
+  }
+
+        if (costRatio > termMaxCostRatio) {
+                //              cout<<"Too big cost ratio "<<costRatio<<endl;
+                return -1;
+        }
+
+#if 0   
+        cout<<
+                "pvs="<<leaf->mPvsSize<<
+                " rays="<<leaf->rays.size()<<
+                " rc="<<leaf->GetAvgRayContribution()<<
+                " axis="<<axis<<endl;
+#endif
+        
+        return axis;
+}
+
+
+                                                        
+
+float
+VspKdTree::EvalCostRatio(
+                                                                                         VspKdTreeLeaf *leaf,
+                                                                                         const int axis,
+                                                                                         const float position,
+                                                                                         int &raysBack,
+                                                                                         int &raysFront,
+                                                                                         int &pvsBack,
+                                                                                         int &pvsFront
+                                                                                         )
+{
+        raysBack = 0;
+        raysFront = 0;
+        pvsFront = 0;
+        pvsBack = 0;
+
+        float newCost;
+
+        Intersectable::NewMail(3);
+        
+        // eval pvs size
+        int pvsSize = leaf->GetPvsSize();
+
+        Intersectable::NewMail(3);
+
+        if (axis <= VspKdTreeNode::SPLIT_Z) {
+    // this is the main ray classification loop!
+    for(VspKdTreeNode::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).ComputeRayIntersection(axis, position, (*ri).mRay->mT);
+                                //                              (*ri).mRay->mSide = side;
+                                
+                                if (side <= 0)
+                                        raysBack++;
+                                
+                                if (side >= 0)
+                                        raysFront++;
+                                
+                                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+      }
+                
+                AxisAlignedBox3 box = GetBBox(leaf);
+                
+                float minBox = box.Min(axis);
+                float maxBox = box.Max(axis);
+                float sizeBox = maxBox - minBox;
+                
+                //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+                float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
+
+                newCost = ct_div_ci + sum/sizeBox;
+                
+  } else {
+                
+                // directional split
+    for(VspKdTreeNode::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);
+
+      }
+                
+                AxisAlignedBox3 box = GetDirBBox(leaf);
+                
+                float minBox = box.Min(axis-3);
+                float maxBox = box.Max(axis-3);
+                float sizeBox = maxBox - minBox;
+                
+                //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+                float sum =
+                        pvsBack*(position - minBox) +
+                        pvsFront*(maxBox - position);
+                //              float sum = pvsBack + pvsFront;
+                newCost = ct_div_ci + sum/sizeBox;
+  }
+
+        //      cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
+        //  float oldCost = leaf->rays.size();
+        float oldCost = pvsSize;
+  float ratio = newCost/oldCost;
+        //      cout<<"ratio="<<ratio<<endl;
+        return ratio;
+}
+
+float
+VspKdTree::BestCostRatioRegular(
+                                                                                                                        VspKdTreeLeaf *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;
+
+        bool onlyDrivingAxis = true; 
+
+        for (axis = 0; axis < 6; axis++) {
+                if (!onlyDrivingAxis || axis == sAxis || axis == dAxis) {
+                        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]
+                                                                                                                                                         );
+                        
+                        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
+VspKdTree::BestCostRatioHeuristic(
+                                                                                                                                VspKdTreeLeaf *leaf,
+                                                                                                                                int &axis,
+                                                                                                                                float &position,
+                                                                                                                                int &raysBack,
+                                                                                                                                int &raysFront,
+                                                                                                                                int &pvsBack,
+                                                                                                                                int &pvsFront
+                                                                                                                                )
+{
+        AxisAlignedBox3 box = GetBBox(leaf);
+        //      AxisAlignedBox3 dirBox = GetDirBBox(node);
+        
+        axis = box.Size().DrivingAxis();
+                
+        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 minBox = box.Min(axis);
+        float maxBox = box.Max(axis);
+        float sizeBox = maxBox - minBox;
+        
+        float minBand = minBox + 0.1*(maxBox - minBox);
+        float maxBand = minBox + 0.9*(maxBox - minBox);
+        
+        float sum = rr*sizeBox;
+        float minSum = 1e20;
+        
+        Intersectable::NewMail();
+        // set all object as belonging to the fron pvs
+        for(VspKdTreeNode::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: {
+                        rl++;
+                        ray = (VssRay *) (*ci).data;
+                        Intersectable *object = ray->mTerminationObject;
+                        if (object && !object->Mailed()) {
+                                object->Mail();
+                                pl++;
+                        }
+                        break;
+                }
+                case SortableEntry::ERayMax: {
+                        rr--;
+                        ray = (VssRay *) (*ci).data;
+                        Intersectable *object = ray->mTerminationObject;
+                        if (object) {
+                                if (--object->mCounter == 0)
+                                        pr--;
+                        }
+                        break;
+                }
+                }
+                if ((*ci).value > minBand && (*ci).value < maxBand) {
+                        
+                        sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);
+                        
+                        //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
+                        //      cout<<"cost= "<<sum<<endl;
+                        
+                        if (sum < minSum) {
+                                minSum = sum;
+                                position = (*ci).value;
+                                
+                                raysBack = rl;
+                                raysFront = rr;
+                                
+                                pvsBack = pl;
+                                pvsFront = pr;
+                                
+      }
     }
-    case SPLIT_SAH: {
-      int objectsBack, objectsFront;
-      float costRatio;
-      bool mOnlyDrivingAxis = false;
-      if (mOnlyDrivingAxis) {
-                                axis = box.Size().DrivingAxis();
-                                costRatio = BestCostRatio(leaf,
-                                                                                                                                        box,
-                                                                                                                                        axis,
-                                                                                                                                        position,
-                                                                                                                                        objectsBack,
-                                                                                                                                        objectsFront);
-      } else {
-                                costRatio = MAX_FLOAT;
-                                for (int i=0; i < 3; i++) {
-                                        float p;
-                                        float r = BestCostRatio(leaf,
-                                                                                                                                        box,
-                                                                                                                                        i,
-                                                                                                                                        p,
-                                                                                                                                        objectsBack,
-                                                                                                                                        objectsFront);
-                                        if (r < costRatio) {
-                                                costRatio = r;
-                                                axis = i;
-                                                position = p;
-                                        }
-                                }
-      }
-      
-      if (costRatio > mMaxCostRatio) {
-                                //      cout<<"Too big cost ratio "<<costRatio<<endl;
-                                axis = -1;
-      }
-      break;
-    }
-    
-    }
-  return axis;
-}
-
-ViewCellKdNode *
-ViewCellKd::SubdivideNode(
-                      ViewCellKdLeaf *leaf,
-                      const AxisAlignedBox3 &box,
-                      AxisAlignedBox3 &backBBox,
-                      AxisAlignedBox3 &frontBBox
-                      )
-{
-  
-  if (TerminationCriteriaMet(leaf))
-    return leaf;
-  
-  float position;
-  
-  // select subdivision axis
-  int axis = SelectPlane( leaf, box, position );
-
-  if (axis == -1) {
-    return leaf;
-  }
-  
-  mStat.nodes+=2;
-  mStat.splits[axis]++;
-  
-  // add the new nodes to the tree
-  ViewCellKdInterior *node = new ViewCellKdInterior(leaf->mParent);
-
-  node->mAxis = axis;
-  node->mPosition = position;
-  node->mBox = box;
-  
-  backBBox = box;
-  frontBBox = box;
-  
-  // first count ray sides
-  int objectsBack = 0;
-  int objectsFront = 0;
-  
-  backBBox.SetMax(axis, position);
-  frontBBox.SetMin(axis, position);
-
-  ObjectContainer::const_iterator mi;
-
-  for ( mi = leaf->mObjects.begin();
-        mi != leaf->mObjects.end();
-        mi++) {
-    // determine the side of this ray with respect to the plane
-    AxisAlignedBox3 box = (*mi)->GetBox();
-    if (box.Max(axis) > position ) 
-      objectsFront++;
-    
-    if (box.Min(axis) < position )
-      objectsBack++;
-  }
-
-  
-  ViewCellKdLeaf *back = new ViewCellKdLeaf(node, objectsBack);
-  ViewCellKdLeaf *front = new ViewCellKdLeaf(node, objectsFront);
-
-  // replace a link from node's parent
-  if (  leaf->mParent )
-    leaf->mParent->ReplaceChildLink(leaf, node);
-
-  // and setup child links
-  node->SetupChildLinks(back, front);
-  
-  for (mi = leaf->mObjects.begin();
-       mi != leaf->mObjects.end();
-       mi++) {
-    // determine the side of this ray with respect to the plane
-    AxisAlignedBox3 box = (*mi)->GetBox();
-
-    if (box.Max(axis) >= position )
-      front->mObjects.push_back(*mi);
-    
-    if (box.Min(axis) < position )
-      back->mObjects.push_back(*mi);
-    
-    mStat.objectRefs -= (int)leaf->mObjects.size();
-    mStat.objectRefs += objectsBack + objectsFront;
-  }
-  
-  delete leaf;
-  return node;
-}
-
-
+  }
+  
+  float oldCost = leaf->GetPvsSize();
+  float newCost = ct_div_ci + minSum/sizeBox;
+  float ratio = newCost/oldCost;
+  
+  //  cout<<"===================="<<endl;
+  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
+  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
+        return ratio;
+}
 
 void
-ViewCellKdTreeStatistics::Print(ostream &app) const
-{
-  app << "===== ViewCellKd statistics ===============\n";
-
-  app << "#N_RAYS Number of rays )\n"
-      << rays <<endl;
-  app << "#N_DOMAINS  ( Number of query domains )\n"
-      << queryDomains <<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_MAXOBJECTREFS  ( Max number of rayRefs / leaf )\n" <<
-    maxObjectRefs << "\n";
-
-  app << "#N_NONEMPTYRAYREFS  ( Number of rayRefs in nonEmpty leaves / non empty leaf )\n" <<
-    rayRefsNonZeroQuery/(double)(Leaves() - zeroQueryNodes) << "\n";
-
-  app << "#N_LEAFDOMAINREFS  ( Number of query domain Refs / leaf )\n" <<
-    objectRefs/(double)Leaves() << "\n";
-
-  //  app << setprecision(4);
-
-  app << "#N_PEMPTYLEAVES  ( Percentage of leaves with zero query domains )\n"<<
-    zeroQueryNodes*100/(double)Leaves()<<endl;
-
-  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
-    maxDepthNodes*100/(double)Leaves()<<endl;
-
-  app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
-    minCostNodes*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 ViewCellKd statistics ==========\n";
-
-}
-
-
-
-void
-ViewCellKd::EvaluateLeafStats(const TraversalData &data)
-{
-
-  // the node became a leaf -> evaluate stats for leafs
-  ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)data.mNode;
-
-  if (data.mDepth > mTermMaxDepth)
-    mStat.maxDepthNodes++;
-  
-  if ( (int)(leaf->mObjects.size()) < mTermMinCost)
-    mStat.minCostNodes++;
-  
-  
-  if ( (int)(leaf->mObjects.size()) > mStat.maxObjectRefs)
-    mStat.maxObjectRefs = (int)leaf->mObjects.size();
-  
-}
-
-
-
-void
-ViewCellKd::SortSplitCandidates(
-                            ViewCellKdLeaf *node,
-                            const int axis
-                            )
-{
+VspKdTree::SortSplitCandidates(
+                                                                                                                 VspKdTreeLeaf *node,
+                                                                                                                 const int axis
+                                                                                                                 )
+{
+  
   splitCandidates->clear();
   
-  int requestedSize = 2 * (int)node->mObjects.size();
+  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>;
@@ -390 +712 @@
   
   splitCandidates->reserve(requestedSize);
-  
+
   // insert all queries 
-  for(ObjectContainer::const_iterator mi = node->mObjects.begin();
-      mi != node->mObjects.end();
-      mi++) {
-    AxisAlignedBox3 box = (*mi)->GetBox();
-
-    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MIN,
-                                             box.Min(axis),
-                                             *mi)
-                               );
+  for(VspKdTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
+      ri < node->rays.end();
+      ri++) {
+    bool positive = (*ri).mRay->HasPosDir(axis);
+    splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
+                                                                                                                                                                                 SortableEntry::ERayMax,
+                                                                                                                                                                                 (*ri).ExtrapOrigin(axis),
+                                                                                                                                                                                 (void *)&*ri)
+                                                                                                                         );
+                
+    splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
+                                                                                                                                                                                 SortableEntry::ERayMin,
+                                                                                                                                                                                 (*ri).ExtrapTermination(axis),
+                                                                                                                                                                                 (void *)&*ri)
+                                                                                                                         );
+  }
+        
+  stable_sort(splitCandidates->begin(), splitCandidates->end());
+}
+
+
+void
+VspKdTree::EvaluateLeafStats(const TraversalData &data)
+{
+
+  // the node became a leaf -> evaluate stats for leafs
+  VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)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 (0 && 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();
+
+}
+
+
+
+VspKdTreeNode *
+VspKdTree::SubdivideNode(
+                                                                                         VspKdTreeLeaf *leaf,
+                                                                                         const AxisAlignedBox3 &box,
+                                                                                         AxisAlignedBox3 &backBBox,
+                                                                                         AxisAlignedBox3 &frontBBox
+                                                                                         )
+{
+  
+  if ( (leaf->GetPvsSize() < termMinPvs) ||
+                         (leaf->rays.size() < termMinRays) ||
+                         //                      (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
+       (leaf->depth >= termMaxDepth) ||
+                         SqrMagnitude(box.Size()) <= termMinSize ) {
+
+#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);
+
+        //      cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
+        //              pvsFront<<endl;
+
+  if (axis == -1) {
+    return leaf;
+  }
+  
+  stat.nodes+=2;
+  stat.splits[axis]++;
+
+  // add the new nodes to the tree
+  VspKdTreeInterior *node = new VspKdTreeInterior(leaf->parent);
+
+  node->axis = axis;
+  node->position = position;
+  node->bbox = box;
+  node->dirBBox = GetDirBBox(leaf);
+  
+  backBBox = box;
+  frontBBox = box;
+
+  VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
+        back->SetPvsSize(pvsBack);
+  VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
+        front->SetPvsSize(pvsFront);
+
+  // 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 <= VspKdTreeNode::SPLIT_Z) {
+                backBBox.SetMax(axis, position);
+    frontBBox.SetMin(axis, position);
+                
+                for(VspKdTreeNode::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();
+
+                                // determine the side of this ray with respect to the plane
+                                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
+
+                                if (side == 0) {
+                                        if ((*ri).mRay->HasPosDir(axis)) {
+                                                back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+                                                                                                                                                                                        (*ri).mMinT,
+                                                                                                                                                                                        (*ri).mRay->mT)
+                                                                                                 );
+                                                front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+                                                                                                                                                                                         (*ri).mRay->mT,
+                                                                                                                                                                                         (*ri).mMaxT));
+                                        } else {
+                                                back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+                                                                                                                                                                                        (*ri).mRay->mT,
+                                                                                                                                                                                        (*ri).mMaxT));
+                                                front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
+                                                                                                                                                                                         (*ri).mMinT,
+                                                                                                                                                                                         (*ri).mRay->mT));
+                                        }
+                                } else
+                                        if (side == 1) 
+                                                front->AddRay(*ri);
+                                        else
+                                                back->AddRay(*ri);
+      } else
+                                (*ri).mRay->Unref();
+    }
+  } else {
+    // rays front/back
     
     
-    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MAX,
-                                             box.Max(axis),
-                                             *mi)
-                               );
-  }
-  
-  stable_sort(splitCandidates->begin(), splitCandidates->end());
-}
-
-
-float
-ViewCellKd::BestCostRatio(
-                                                                                        ViewCellKdLeaf *node,
-                                                                                        const AxisAlignedBox3 &box,
-                                                                                        const int axis,
-                                                                                        float &position,
-                                                                                        int &objectsBack,
-                                                                                        int &objectsFront
-                                                                                        )
-{
-
-  SortSplitCandidates(node, axis);
-  
-  // go through the lists, count the number of objects left and right
-  // and evaluate the following cost funcion:
-  // C = ct_div_ci  + (ol + or)/queries
-
-  float totalIntersections = 0.0f;
-  vector<SortableEntry>::const_iterator ci;
-
-  for(ci = splitCandidates->begin();
-      ci < splitCandidates->end();
-      ci++) 
-    if ((*ci).type == SortableEntry::BOX_MIN) {
-      totalIntersections += (*ci).intersectable->IntersectionComplexity();
+    for(VspKdTreeNode::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();
     }
-        
-  float intersectionsLeft = 0;
-  float intersectionsRight = totalIntersections;
-        
-  int objectsLeft = 0, objectsRight = (int)node->mObjects.size();
-  
-  float minBox = box.Min(axis);
-  float maxBox = box.Max(axis);
-  float boxArea = box.SurfaceArea();
-  
-  float minBand = minBox + mSplitBorder*(maxBox - minBox);
-  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
-  
-  float minSum = 1e20;
-  
-  for(ci = splitCandidates->begin();
-      ci < splitCandidates->end();
-      ci++) {
-    switch ((*ci).type) {
-    case SortableEntry::BOX_MIN:
-      objectsLeft++;
-      intersectionsLeft += (*ci).intersectable->IntersectionComplexity();
-      break;
-    case SortableEntry::BOX_MAX:
-      objectsRight--;
-      intersectionsRight -= (*ci).intersectable->IntersectionComplexity();
-      break;
-    }
-
-    if ((*ci).value > minBand && (*ci).value < maxBand) {
-      AxisAlignedBox3 lbox = box;
-      AxisAlignedBox3 rbox = box;
-      lbox.SetMax(axis, (*ci).value);
-      rbox.SetMin(axis, (*ci).value);
-
-      float sum;
-      if (mSahUseFaces)
-                                sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
-      else
-                                sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
+  }
+        
+  // update stats
+  stat.rayRefs -= leaf->rays.size();
+  stat.rayRefs += raysBack + raysFront;
+
+  
+  delete leaf;
+  return node;
+}
+
+
+
+
+
+
+int
+VspKdTree::ReleaseMemory(const int time)
+{
+  stack<VspKdTreeNode *> tstack;
+  
+  // find a node in the tree which subtree will be collapsed
+  int maxAccessTime = time - accessTimeThreshold;
+  int released;
+
+  tstack.push(root);
+
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
+    if (!node->IsLeaf()) {
+      VspKdTreeInterior *in = (VspKdTreeInterior *)node;
+      //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
+      if (in->depth >= minCollapseDepth &&
+          in->lastAccessTime <= maxAccessTime) {
+        released = CollapseSubtree(node, time);
+        break;
+      }
       
-      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
-      //      cout<<"cost= "<<sum<<endl;
-      
-      if (sum < minSum) {
-                                minSum = sum;
-                                position = (*ci).value;
-                                
-                                objectsBack = objectsLeft;
-                                objectsFront = objectsRight;
+      if (in->back->GetAccessTime() <  
+          in->front->GetAccessTime()) {
+        tstack.push(in->front);
+        tstack.push(in->back);
+      } else {
+        tstack.push(in->back);
+        tstack.push(in->front);
       }
     }
   }
-  
-  float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
-  float newCost = mCt_div_ci + minSum/boxArea;
-  float ratio = newCost/oldCost;
-  
-#if 0
-  cout<<"===================="<<endl;
-  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
-      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
+
+  while (tstack.empty()) {
+    // could find node to collaps...
+    //    cout<<"Could not find a node to release "<<endl;
+    break;
+  }
+  
+  return released;
+}
+
+
+
+
+VspKdTreeNode *
+VspKdTree::SubdivideLeaf(
+                                                                                         VspKdTreeLeaf *leaf,
+                                                                                         const float sizeThreshold
+                                                                                         )
+{
+  VspKdTreeNode *node = leaf;
+
+  AxisAlignedBox3 leafBBox = GetBBox(leaf);
+
+        static int pass = 0;
+        pass ++;
+        
+  // check if we should perform a dynamic subdivision of the leaf
+  if (
+                        //      leaf->rays.size() > (unsigned)termMinCost &&
+                        (leaf->GetPvsSize() >= termMinPvs) &&
+      SqrMagnitude(leafBBox.Size()) > sizeThreshold) {
+    
+                // memory check and realese...
+    if (GetMemUsage() > maxTotalMemory) {
+      ReleaseMemory( pass );
+    }
+    
+    AxisAlignedBox3 backBBox, frontBBox;
+
+    // subdivide the node
+    node =
+      SubdivideNode(leaf,
+                                                                                leafBBox,
+                                                                                backBBox,
+                                                                                frontBBox
+                                                                                );
+  }
+  return node;
+}
+
+
+
+void
+VspKdTree::UpdateRays(VssRayContainer &remove,
+                                                                                VssRayContainer &add
+                                                                                )
+{
+  VspKdTreeLeaf::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++) {
+    AddRay(*ri);
+  }
+  
+
+}
+
+
+void
+VspKdTree::RemoveRay(VssRay *ray,
+                                                                         vector<VspKdTreeLeaf *> *affectedLeaves,
+                                                                         const bool removeAllScheduledRays
+                                                                         )
+{
+        
+  stack<RayTraversalData> tstack;
+
+  tstack.push(RayTraversalData(root, VspKdTreeNode::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...
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)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
+VspKdTree::AddRay(VssRay *ray)
+{
+
+  stack<RayTraversalData> tstack;
+  
+  tstack.push(RayTraversalData(root, VspKdTreeNode::RayInfo(ray)));
+  
+  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...
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)data.node;
+      leaf->AddRay(data.rayData);
+      stat.addedRayRefs++;
+    }
+  }
+}
+
+void
+VspKdTree::TraverseInternalNode(
+                                                                                                                        RayTraversalData &data,
+                                                                                                                        stack<RayTraversalData> &tstack)
+{
+  VspKdTreeInterior *in =  (VspKdTreeInterior *) data.node;
+  
+  if (in->axis <= VspKdTreeNode::SPLIT_Z) {
+
+    // determine the side of this ray with respect to the plane
+    int side = in->ComputeRayIntersection(data.rayData,
+                                                                                                                                                                        data.rayData.mRay->mT);
+    
+  
+    if (side == 0) {
+      if (data.rayData.mRay->HasPosDir(in->axis)) {
+                                tstack.push(RayTraversalData(in->back,
+                                                                                                                                                 VspKdTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                                                                                                                                        data.rayData.mMinT,
+                                                                                                                                                                                                                                        data.rayData.mRay->mT))
+                                                                                );
+                                
+                                tstack.push(RayTraversalData(in->front,
+                                                                                                                                                 VspKdTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                                                                                                                                        data.rayData.mRay->mT,
+                                                                                                                                                                                                                                        data.rayData.mMaxT
+                                                                                                                                                                                                                                        ))
+                                                                                );
+        
+      } else {
+                                tstack.push(RayTraversalData(in->back,
+                                                                                                                                                 VspKdTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                                                                                                                                        data.rayData.mRay->mT,
+                                                                                                                                                                                                                                        data.rayData.mMaxT
+                                                                                                                                                                                                                                        ))
+                                                                                );
+                                
+                                tstack.push(RayTraversalData(in->front,
+                                                                                                                                                 VspKdTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                                                                                                                                        data.rayData.mMinT,
+                                                                                                                                                                                                                                        data.rayData.mRay->mT))
+                                                                                );
+                                
+                                
+      }
+    } 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->front, data.rayData));
+                else
+                        tstack.push(RayTraversalData(in->back, data.rayData));
+  }
+}
+
+
+int
+VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
+{
+  // first count all rays in the subtree
+  // use mail 1 for this purpose
+  stack<VspKdTreeNode *> 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
-  return ratio;
-}
-
-int
-ViewCellKd::CastRay(
-                                                                Ray &ray
-                                                                )
-{
-  int hits = 0;
-  
-  stack<RayTraversalData> tStack;
-  
-  float maxt = 1e6;
-  float mint = 0;
-
-  Intersectable::NewMail();
-
-  if (!mBox.GetMinMaxT(ray, &mint, &maxt))
-    return 0;
-
-  if (mint < 0)
-    mint = 0;
-  
-  maxt += Limits::Threshold;
-  
-  Vector3 entp = ray.Extrap(mint);
-  Vector3 extp = ray.Extrap(maxt);
-  
-  ViewCellKdNode *node = mRoot;
-  ViewCellKdNode *farChild;
-  float position;
-  int axis;
-  
-  while (1) {
-    if (!node->IsLeaf()) {
-      ViewCellKdInterior *in = (ViewCellKdInterior *) node;
-      position = in->mPosition;
-      axis = in->mAxis;
-
-      if (entp[axis] <= position) {
-                                if (extp[axis] <= position) {
-                                        node = in->mBack;
-                                        // cases N1,N2,N3,P5,Z2,Z3
-                                        continue;
-                                } else {
-                                        // case N4
-                                        node = in->mBack;
-                                        farChild = in->mFront;
-                                }
-      }
-      else {
-                                if (position <= extp[axis]) {
-                                        node = in->mFront;
-                                        // cases P1,P2,P3,N5,Z1
-                                        continue;
-                                } else {
-                                        node = in->mFront;
-                                        farChild = in->mBack;
-                                        // case P4
-                                }
-                        }
-      // $$ modification 3.5.2004 - hints from Kamil Ghais
-      // case N4 or P4
-      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
-      tStack.push(RayTraversalData(farChild, extp, maxt));
-      extp = ray.GetLoc() + ray.GetDir()*tdist;
-      maxt = tdist;
-                } else {
-                        // compute intersection with all objects in this leaf
-                        ViewCellKdLeaf *leaf = (ViewCellKdLeaf *) node;
-//                      if (ray.mFlags & Ray::STORE_KDLEAVES)
-//                              ray.kdLeaves.push_back(leaf);
-                        
-                        ObjectContainer::const_iterator mi;
-                        for ( mi = leaf->mObjects.begin();
-                                                mi != leaf->mObjects.end();
-                                                mi++) {
-                                Intersectable *object = *mi;
-                                if (!object->Mailed() ) {
-                                        object->Mail();
-                                        if (ray.mFlags & Ray::STORE_TESTED_OBJECTS)
-                                                ray.testedObjects.push_back(object);
-                                        hits += object->CastRay(ray);
-                                }
-                        }
-                        
-                        if (hits && ray.GetType() == Ray::LOCAL_RAY)
-                                if (ray.intersections[0].mT <= maxt)
-                                        break;
-                        
-                        // get the next node from the stack
-                        if (tStack.empty())
-                                break;
-                        
-                        entp = extp;
-                        mint = maxt;
-                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
-                                break;
-                        
-                        RayTraversalData &s  = tStack.top();
-                        node = s.mNode;
-                        extp = s.mExitPoint;
-                        maxt = s.mMaxT;
-                        tStack.pop();
-                }
-  }
-  return hits;
-}
-
-void
-ViewCellKd::CollectObjects(ViewCellKdNode *n, ObjectContainer &objects)
-{
-  stack<ViewCellKdNode *> nodeStack;
-
-  nodeStack.push(n);
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
+
+//    tstat.collapsedSubtrees++;
+//    tstat.collapseDepths += (int)sroot->depth;
+//    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
+  
+  tstack.push(sroot);
+  VssRay::NewMail();
+        
+  while (!tstack.empty()) {
+    collapsedNodes++;
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+
     if (node->IsLeaf()) {
-      ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)node;
-      for (int j=0; j < leaf->mObjects.size(); j++) {
-                                Intersectable *object = leaf->mObjects[j];
-                                if (!object->Mailed()) {
-                                        object->Mail();
-                                        objects.push_back(object);
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
+      for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                                        ri != leaf->rays.end();
+                                        ri++) {
+                                
+                                totalRayCount++;
+                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
+                                        (*ri).mRay->Mail();
+                                        rayCount++;
                                 }
       }
     } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      nodeStack.push(interior->mFront);
-      nodeStack.push(interior->mBack);
+      tstack.push(((VspKdTreeInterior *)node)->back);
+      tstack.push(((VspKdTreeInterior *)node)->front);
     }
   }
-}
-
-// Find random neighbor which was not mailed
-ViewCellKdNode *
-ViewCellKd::FindRandomNeighbor(ViewCellKdNode *n,
-                                                                                                         bool onlyUnmailed
-                                                                                                         )
-{
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-
-  AxisAlignedBox3 box = GetBox(n);
-  int mask = rand();
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
+  
+  VssRay::NewMail();
+
+  // create a new node that will hold the rays
+  VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf( sroot->parent, rayCount );
+  if (  newLeaf->parent )
+    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
+  
+
+  tstack.push( sroot );
+
+  while (!tstack.empty()) {
+
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+
     if (node->IsLeaf()) {
-      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
-        return node;
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      if (interior->mPosition > box.Max(interior->mAxis))
-        nodeStack.push(interior->mBack);
-      else
-        if (interior->mPosition < box.Min(interior->mAxis))
-          nodeStack.push(interior->mFront);
-        else {
-          // random decision
-          if (mask&1)
-            nodeStack.push(interior->mBack);
-          else
-            nodeStack.push(interior->mFront);
-          mask = mask>>1;
-        }
-    }
-  }
-  
-  return NULL;
-}
-
-int
-ViewCellKd::FindNeighbors(ViewCellKdNode *n,
-                      vector<ViewCellKdNode *> &neighbors,
-                      bool onlyUnmailed
-                      )
-{
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-
-  AxisAlignedBox3 box = GetBox(n);
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      if ( node != n && (!onlyUnmailed || !node->Mailed()) ) 
-        neighbors.push_back(node);
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      if (interior->mPosition > box.Max(interior->mAxis))
-                                nodeStack.push(interior->mBack);
-      else
-                                if (interior->mPosition < box.Min(interior->mAxis))
-                                        nodeStack.push(interior->mFront);
-                                else {
-                                        // random decision
-                                        nodeStack.push(interior->mBack);
-                                        nodeStack.push(interior->mFront);
-                                }
-    }
-  }
-  
-  return (int)neighbors.size();
-}
-
-// Find random neighbor which was not mailed
-ViewCellKdNode *
-ViewCellKd::GetRandomLeaf(const Plane3 &plane)
-{
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-  
-  int mask = rand();
-  
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      return node;
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      ViewCellKdNode *next;
-        if (GetBox(interior->mBack).Side(plane) < 0)
-          next = interior->mFront;
-        else
-          if (GetBox(interior->mFront).Side(plane) < 0)
-            next = interior->mBack;
-          else {
-            // random decision
-            if (mask&1)
-              next = interior->mBack;
-            else
-              next = interior->mFront;
-            mask = mask>>1;
-          }
-        nodeStack.push(next);
-    }
-  }
-  
-  
-  return NULL;
-}
-
-void
-ViewCellKd::CollectLeaves(vector<ViewCellKdLeaf *> &leaves)
-{
-  stack<ViewCellKdNode *> nodeStack;
-  nodeStack.push(mRoot);
-
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)node;
-      leaves.push_back(leaf);
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      nodeStack.push(interior->mBack);
-      nodeStack.push(interior->mFront);
-    }
-  }
-}
-
-
-int
-ViewCellKd::CollectLeafPvs()
-{
-  int totalPvsSize = 0;
-  stack<ViewCellKdNode *> nodeStack;
-  
-  nodeStack.push(mRoot);
-  
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
-    if (node->IsLeaf()) {
-      ViewCellKdLeaf *leaf = (ViewCellKdLeaf *)node;
-      for (int j=0; j < leaf->mObjects.size(); j++) {
-        Intersectable *object = leaf->mObjects[j];
-        if (!object->Mailed()) {
-          object->Mail();
-          // add this node to pvs of all nodes it can see
-         /* ViewCellKdPvsMap::iterator ni = object->mKdPvs.mEntries.begin();
-          for (; ni != object->mKdPvs.mEntries.end(); ni++) {
-            ViewCellKdNode *node = (*ni).first;
-            // $$ JB TEMPORARY solution -> should add object PVS or explictly computed
-            // kd tree PVS
-            if (leaf->mKdPvs.AddSample(node))
-              totalPvsSize++;
-          }*/
-        }
+      VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
+      
+      for(VspKdTreeNode::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 {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-      nodeStack.push(interior->mFront);
-      nodeStack.push(interior->mBack);
+      tstack.push(((VspKdTreeInterior *)node)->back);
+      tstack.push(((VspKdTreeInterior *)node)->front);
     }
   }
-
-  return totalPvsSize;
-}
-
-
-ViewCellKdNode *
-ViewCellKd::GetRandomLeaf(const bool onlyUnmailed)
-{
-        stack<ViewCellKdNode *> nodeStack;
-  nodeStack.push(mRoot);
-
-        int mask = rand();
-        
-  while (!nodeStack.empty()) {
-    ViewCellKdNode *node = nodeStack.top();
-    nodeStack.pop();
+  
+  // delete the node and all its children
+  delete sroot;
+  
+//   for(VspKdTreeNode::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
+VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
+{
+        stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+
+        Intersectable::NewMail();
+        int pvsSize = 0;
+        
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+    
+ 
     if (node->IsLeaf()) {
-      if ( (!onlyUnmailed || !node->Mailed()) )
-                                return node;
-    } else {
-      ViewCellKdInterior *interior = (ViewCellKdInterior *)node;
-                        // random decision
-                        if (mask&1)
-                                nodeStack.push(interior->mBack);
-                        else
-                                nodeStack.push(interior->mFront);
-                        mask = mask>>1;
+                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+                        for(VspKdTreeNode::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 {
+                        VspKdTreeInterior *in = (VspKdTreeInterior *)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 NULL;
-}
+        return pvsSize;
+}
+
+void
+VspKdTree::GetRayContributionStatistics(
+                                                                                                                                                        float &minRayContribution,
+                                                                                                                                                        float &maxRayContribution,
+                                                                                                                                                        float &avgRayContribution
+                                                                                                                                                        )
+{
+        stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+        
+        minRayContribution = 1.0f;
+        maxRayContribution = 0.0f;
+        float sumRayContribution = 0.0f;
+        int leaves = 0;
+
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+                
+    if (node->IsLeaf()) {
+                        leaves++;
+                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+                        float c = leaf->GetAvgRayContribution();
+                        if (c > maxRayContribution)
+                                maxRayContribution = c;
+                        if (c < minRayContribution)
+                                minRayContribution = c;
+                        sumRayContribution += c;
+                        
+                } else {
+                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
+                        // both nodes for directional splits
+                        tstack.push(in->front);
+                        tstack.push(in->back);
+                }
+        }
+        
+        cout<<"sum="<<sumRayContribution<<endl;
+        cout<<"leaves="<<leaves<<endl;
+        avgRayContribution = sumRayContribution/(float)leaves;
+}
+
+
+int
+VspKdTree::GenerateRays(const float ratioPerLeaf,
+                                                                                        SimpleRayContainer &rays)
+{
+        stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+        
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+                
+    if (node->IsLeaf()) {
+                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+                        float c = leaf->GetAvgRayContribution();
+                        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 = 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));
+                        }
+                        
+                } else {
+                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
+                        // both nodes for directional splits
+                        tstack.push(in->front);
+                        tstack.push(in->back);
+                }
+        }
+
+        return rays.size();
+}
+
+
+float
+VspKdTree::GetAvgPvsSize()
+{
+        stack<VspKdTreeNode *> tstack;
+  tstack.push(root);
+
+        int sumPvs = 0;
+        int leaves = 0;
+  while (!tstack.empty()) {
+    VspKdTreeNode *node = tstack.top();
+    tstack.pop();
+                
+    if (node->IsLeaf()) {
+                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
+                        // update pvs size
+                        leaf->UpdatePvsSize();
+                        sumPvs += leaf->GetPvsSize();
+                        leaves++;
+                } else {
+                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
+                        // both nodes for directional splits
+                        tstack.push(in->front);
+                        tstack.push(in->back);
+                }
+        }
+
+
+        return sumPvs/(float)leaves;
+}

trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.h

@@ -1 @@
-#ifndef _ViewCellKdTree_H__
-#define _ViewCellKdTree_H__
-
+// ================================================================
+// $Id$
+// ****************************************************************
+// 
+// Initial coding by
+/**
+   @author Jiri Bittner
+*/
+
+#ifndef __VSPKDTREE_H__
+#define __VSPKDTREE_H__
+
+// Standard headers
+#include <iomanip>
+#include <vector>
 #include <functional>
-using namespace std;
-
-#include "Containers.h"
+#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 "Pvs.h"
-
-  
-class ViewCellKdNode;
-class ViewCellKdLeaf;
-class ViewCellKdInterior;
-class Intersectable;
-
 
 // --------------------------------------------------------------
 // Static statistics for kd-tree search
 // --------------------------------------------------------------
-class ViewCellKdTreeStatistics 
+class VspKdStatistics:  public StatisticsBase
 {
 public:  
@@ -29 +43 @@
   // totals number of rays
   int rays;
-  // total number of query domains
+        // initial size of the pvs
+  int initialPvsSize;
+        // total number of query domains
   int queryDomains;
   // total number of ray references
   int rayRefs;
-  // refs in non empty leafs
-  int rayRefsNonZeroQuery;
-  // total number of query references
-  int objectRefs;
-  // nodes with zero queries
-  int zeroQueryNodes;
-  // max depth nodes
+
+        // max depth nodes
   int maxDepthNodes;
   // max depth nodes
-  int minCostNodes;
+  int minPvsNodes;
+  int minRaysNodes;
+        // max ray contribution nodes
+  int maxRayContribNodes;
+        // max depth nodes
+  int minSizeNodes;
+        
   // max number of rays per node
-  int maxObjectRefs;
+  int maxRayRefs;
   // number of dynamically added ray refs
   int addedRayRefs;
@@ -51 +68 @@
   
   // Constructor
-  ViewCellKdTreeStatistics() {
+  VspKdStatistics() {
     Reset();
   }
@@ -64 +81 @@
       splits[i] = 0;
     rays = queryDomains = 0;
-    rayRefs = rayRefsNonZeroQuery = objectRefs = 0;
-    zeroQueryNodes = 0;
+    rayRefs = 0;
     maxDepthNodes = 0;
-    minCostNodes = 0;
-    maxObjectRefs = 0;
+    minPvsNodes = 0;
+    minRaysNodes = 0;
+    maxRayRefs = 0;
     addedRayRefs = removedRayRefs = 0;
-  }
-
+                initialPvsSize = 0;
+                maxRayContribNodes = 0;
+                minSizeNodes = 0;
+  }
+
+  
   void
   Print(ostream &app) const;
 
-  friend ostream &operator<<(ostream &s, const ViewCellKdTreeStatistics &stat) {
+  friend ostream &operator<<(ostream &s, const VspKdStatistics &stat) {
     stat.Print(s);
     return s;
@@ -82 +103 @@
 };
 
-  
-class ViewCellKdInterior;
-/** Abstract class for kd-tree node */
-class ViewCellKdNode {
+
+// --------------------------------------------------------------
+// For sorting rays
+// --------------------------------------------------------------
+struct  SortableEntry
+{
+  enum EType {
+    ERayMin,
+    ERayMax
+  };
+
+  int type;
+  float value;
+  void *data;
+  
+  SortableEntry() {}
+  SortableEntry(const int t, const float v, void *d):type(t),
+                                                                                                                                                                                                                 value(v),
+                                                                                                                                                                                                                 data(d) {}
+        
+  friend bool operator<(const SortableEntry &a, const SortableEntry &b) {
+    return a.value < b.value;
+  }
+};
+
+
+class VspKdTreeInterior;
+
+
+// --------------------------------------------------------------
+// KD-tree node - abstract class
+// --------------------------------------------------------------
+class VspKdTreeNode
+{
+public:
+
+#define USE_FIXEDPOINT_T 1
+
+struct RayInfo
+{
+        // pointer to the actual ray
+        VssRay *mRay;
+        
+        // endpoints  - do we need them?
+#if USE_FIXEDPOINT_T
+        short mMinT, mMaxT;
+#else
+        float mMinT, mMaxT;
+#endif
+        
+        RayInfo():mRay(NULL) {}
+        
+        RayInfo(VssRay *r):mRay(r), mMinT(0),
+#if USE_FIXEDPOINT_T
+#define FIXEDPOINT_ONE 0x7FFE
+                                                                                                //                        mMaxT(0xFFFF)
+                                                                                                                                                        mMaxT(FIXEDPOINT_ONE)
+#else
+                mMaxT(1.0f)
+#endif
+        {}
+        
+        RayInfo(VssRay *r,
+                                                                                const float _min,
+                                                                                const float _max
+                                                                                ):mRay(r) {
+                SetMinT(_min);
+                SetMaxT(_max);
+        }
+        
+        RayInfo(VssRay *r,
+                                                                                const short _min,
+                                                                                const float _max
+                                                                                ):mRay(r), mMinT(_min) {
+                SetMaxT(_max);
+        }
+        
+        RayInfo(VssRay *r,
+                                                                                const float _min,
+                                                                                const short _max
+                                                                                ):mRay(r), mMaxT(_max) {
+                SetMinT(_min);
+        }
+        
+        friend bool operator<(const RayInfo &a, const RayInfo &b) {
+                return a.mRay < b.mRay;
+        }
+  
+        
+        float ExtrapOrigin(const int axis) const {
+                return mRay->GetOrigin(axis) + GetMinT()*mRay->GetDir(axis);
+        }
+        
+        float ExtrapTermination(const int axis) const {
+                return mRay->GetOrigin(axis) + GetMaxT()*mRay->GetDir(axis);
+        }
+        
+#if USE_FIXEDPOINT_T
+        float GetMinT () const { return mMinT/(float)(FIXEDPOINT_ONE); }
+        float GetMaxT () const { return mMaxT/(float)(FIXEDPOINT_ONE); }
+        
+        void SetMinT (const float t) {
+                mMinT = (short) (t*(float)(FIXEDPOINT_ONE));
+        }
+        
+        void SetMaxT (const float t) {
+                mMaxT = (short) (t*(float)(FIXEDPOINT_ONE));
+                mMaxT++;
+                //      if (mMaxT!=0xFFFF)
+                //      mMaxT++;
+        }
+#else
+        float GetMinT () const { return mMinT; }
+        float GetMaxT () const { return mMaxT; }
+        
+        void SetMinT (const float t) { mMinT = t; }
+        void SetMaxT (const float t) { mMaxT = t; }
+#endif 
+
+
+  int ComputeRayIntersection(const float axis,
+                                                                                                                 const float position,
+                                                                                                                 float &t
+                                                                                                                 ) const {
+                
+                // intersect the ray with the plane
+    float denom = mRay->GetDir(axis);
+    
+    if (fabs(denom) < 1e-20)
+      //if (denom == 0.0f)
+      return (mRay->GetOrigin(axis) > position) ? 1 : -1;
+    
+    t = (position - mRay->GetOrigin(axis))/denom;
+
+    if (t < GetMinT())
+      return (denom > 0) ? 1 : -1;
+
+    if (t > GetMaxT())
+      return (denom > 0) ? -1 : 1;
+
+                return 0;
+        }
+
+
+};
+
+
+        typedef vector<RayInfo> RayInfoContainer;
+        
+  enum { EInterior, ELeaf };
+
+  /////////////////////////////////
+  // The actual data goes here
+  
+  // link to the parent
+  VspKdTreeInterior *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 VspKdTreeNode(VspKdTreeInterior *p);
+
+  
+  virtual ~VspKdTreeNode() {};
+  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 VspKdTreeInterior :
+  public VspKdTreeNode
+{
+public:
+  // plane in local modelling coordinates
+  float position;
+
+  // pointers to children
+  VspKdTreeNode *back, *front;
+
+  // the bbox of the node
+  AxisAlignedBox3 bbox;
+
+  // the bbox of the node
+  AxisAlignedBox3 dirBBox;
+  
+  // data for caching
+  long accesses;
+  long lastAccessTime;
+  
+  VspKdTreeInterior(VspKdTreeInterior *p):VspKdTreeNode(p),
+                                        back(NULL),
+                                        front(NULL),
+                                        accesses(0),
+                                        lastAccessTime(-1)
+  { }
+
+  virtual int GetAccessTime() {
+    return lastAccessTime;
+  }
+
+  void SetupChildLinks(VspKdTreeNode *b, VspKdTreeNode *f) {
+    back = b;
+    front = f;
+    b->parent = f->parent = this;
+  }
+
+  void ReplaceChildLink(VspKdTreeNode *oldChild, VspKdTreeNode *newChild) {
+    if (back == oldChild)
+      back = newChild;
+    else
+      front = newChild;
+  }
+
+  virtual int Type() const  { return EInterior; }
+  
+  virtual ~VspKdTreeInterior() {
+    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 ComputeRayIntersection(const RayInfo &rayData,
+                                                                                                                 float &t
+                                                                                                                 ) {
+                return rayData.ComputeRayIntersection(axis, position, t);
+  }
+
+};
+
+
+// --------------------------------------------------------------
+// KD-tree node - leaf node
+// --------------------------------------------------------------
+class VspKdTreeLeaf :
+  public VspKdTreeNode
+{
+private:
+        int mPvsSize;
 public:
   static int mailID;
   int mailbox;
   
+  RayInfoContainer rays;
+        bool mValidPvs;
+        
+  VspKdTreeLeaf(VspKdTreeInterior *p,
+                                                        const int nRays
+                                                        ):VspKdTreeNode(p), rays(), mPvsSize(0), mValidPvs(false) {
+    rays.reserve(nRays);
+  }
+  
+  virtual ~VspKdTreeLeaf() { }
+
+  virtual int Type() const  { return ELeaf; }
+
+  virtual void Print(ostream &s) const {
+    s<<endl<<"L: r="<<rays.size()<<endl;
+  };
+  
+  void AddRay(const RayInfo &data) {
+                mValidPvs = false;
+    rays.push_back(data);
+    data.mRay->Ref();
+  }
+        
+        int GetPvsSize() const {
+                return mPvsSize;
+        }
+        void SetPvsSize(const int s) {
+                mPvsSize = s;
+        }
+
+        void
+        UpdatePvsSize();
+
   void Mail() { mailbox = mailID; }
   static void NewMail() { mailID++; }
   bool Mailed() const { return mailbox == mailID; }
-
-
-  ViewCellKdNode(ViewCellKdInterior *parent);
-
-  /** Determines whether this node is a leaf or interior node
-      @return true if leaf
-  */
-  virtual bool IsLeaf() const = 0;
-
-  /** Determines whether this node is the root of the tree
-      @return true if root
-  */
-  virtual bool IsRoot() const {
-    return mParent == NULL;
-  }
-
-  /** Parent of the node - the parent is a little overhead for maintanance of the tree,
-      but allows various optimizations of tree traversal algorithms */
-  ViewCellKdInterior *mParent;
-  int mDepth;
+  
+  bool Mailed(const int mail) {
+    return mailbox >= mailID + mail;
+  }
+
+        float GetAvgRayContribution() const {
+                return GetPvsSize()/((float)rays.size() + Limits::Small);
+        }
+
+        float GetSqrRayContribution() const {
+                return sqr(GetPvsSize()/((float)rays.size() + Limits::Small));
+        }
+
 };
 
-/** Implementation of the kd-tree interior node */
-class ViewCellKdInterior : public ViewCellKdNode {
-
-public:
-    
-  ViewCellKdInterior(ViewCellKdInterior *parent):ViewCellKdNode(parent), mBack(NULL), mFront(NULL) {}
-
-  /** \sa ViewCellKdNode::IsLeaf() */
-  virtual bool IsLeaf() const { return false; }
-
-  /** splitting axis */
-  int mAxis;
-  /** splitting position, absolute position within the bounding box of this node */
-  float mPosition;
-  /** bounding box of interior node */
-  AxisAlignedBox3 mBox;
-
-  /** back node */
-  ViewCellKdNode *mBack;
-  /** front node */
-  ViewCellKdNode *mFront;
-
-  void SetupChildLinks(ViewCellKdNode *b, ViewCellKdNode *f) {
-    mBack = b;
-    mFront = f;
-    b->mParent = f->mParent = this;
-  }
-  
-  void ReplaceChildLink(ViewCellKdNode *oldChild, ViewCellKdNode *newChild) {
-    if (mBack == oldChild)
-      mBack = newChild;
-    else
-      mFront = newChild;
-  }
-  
-  
-};
-  
-  
-/** Implementation of the kd-tree leaf node */
-class ViewCellKdLeaf : public ViewCellKdNode {
-public:
-  ViewCellKdLeaf(ViewCellKdInterior *parent, const int objects):ViewCellKdNode(parent) {
-    mObjects.reserve(objects);
-  }
-  
-  void AddPassingRay(const Ray &ray,
-                                                                                 const int contributions) {
-    mPassingRays.AddRay(ray, contributions);
-                //              Debug << "adding passing ray" << endl;
-  }
-
-        
-        void AddPassingRay2(const Ray &ray,
-                                                                                        const int objects,
-                                                                                        const int viewcells
-                                                                                        ) {
-    mPassingRays.AddRay2(ray, objects, viewcells);
-                //              Debug << "adding passing ray" << endl;
-  }
-
-  /** \sa ViewCellKdNode::IsLeaf() */
-  virtual bool IsLeaf() const { return true; }
-  
-
-  
-
-  /** pointers to occluders contained in this node */
-  ObjectContainer mObjects;
-  
-  /** pointers to viewcells contained in this node */
-  //  ViewCellContainer mViewCells;
-
-  /** Ray set description of the rays passing through this node */
-  PassingRaySet mPassingRays;
-        
-  /** PVS consisting of visible ViewCellKd nodes */
-  ViewCellKdPvs mKdPvs;
-        
-        /** PVS consisting of visible objects */
-  ViewCellPvs mPvs;
-};
-
-  
-
-/** ViewCellKd for indexing scene entities - occluders/occludees/viewcells */
-class ViewCellKd {
-
-protected:
+// Inline functions
+inline
+VspKdTreeNode::VspKdTreeNode(VspKdTreeInterior *p):
+  parent(p), axis(-1), depth(p ? p->depth + 1 : 0) {}
+
+
+
+// ---------------------------------------------------------------
+// Main LSDS search class
+// ---------------------------------------------------------------
+class VspKdTree
+{
   struct TraversalData
   {  
-    ViewCellKdNode *mNode;
-    AxisAlignedBox3 mBox;
-    int mDepth;
-    float mPriority;
+    VspKdTreeNode *node;
+    AxisAlignedBox3 bbox;
+    int depth;
+    float priority;
     
     TraversalData() {}
 
-    TraversalData(ViewCellKdNode *n, const float p):
-      mNode(n), mPriority(p)
+    TraversalData(VspKdTreeNode *n, const float p):
+      node(n), priority(p)
     {}
-    
-    TraversalData(ViewCellKdNode *n,
+
+    TraversalData(VspKdTreeNode *n,
                    const AxisAlignedBox3 &b,
                    const int d):
-      mNode(n), mBox(b), mDepth(d) {}
+      node(n), bbox(b), depth(d) {}
     
-
-    bool operator<(
-                   const TraversalData &b) const {
-      ViewCellKdLeaf *leafa = (ViewCellKdLeaf *) mNode;
-      ViewCellKdLeaf *leafb = (ViewCellKdLeaf *) b.mNode;
-      return 
-        leafa->mObjects.size()*mBox.SurfaceArea()
-        <
-        leafb->mObjects.size()*b.mBox.SurfaceArea();
-    }
-
-
+                
     // comparator for the 
     struct less_priority : public
     binary_function<const TraversalData, const TraversalData, bool> {
-      
+                        
       bool operator()(const TraversalData a, const TraversalData b) {
-                                return a.mPriority < b.mPriority;
+                                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();
+      VspKdTreeLeaf *leafa = (VspKdTreeLeaf *) a.node;
+      VspKdTreeLeaf *leafb = (VspKdTreeLeaf *) b.node;
+#if 0
+                        return
+                                leafa->rays.size()*a.bbox.GetVolume()
+                                <
+                                leafb->rays.size()*b.bbox.GetVolume();
+#endif
+#if 1
+                        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 0
+                        return
+                                leafa->GetPvsSize()*leafa->rays.size()
+                                <
+                                leafb->GetPvsSize()*leafb->rays.size();
+#endif
+    }
   };
-
-  
-
-public:
-
-  enum {SPLIT_OBJECT_MEDIAN,
-        SPLIT_SPATIAL_MEDIAN,
-        SPLIT_SAH};
-  
-  ViewCellKd();
-  
+  
+  // simplified data for ray traversal only...
+
+  struct RayTraversalData {
     
-  /** Insert view cell into the tree */
-  virtual void InsertViewCell(ViewCell *viewCell) {
-    //      mRoot->mViewcells.push_back(viewCell);
-  }
-
-  virtual bool Construct();
-
-  /** Check whether subdivision criteria are met for the given subtree.
-      If not subdivide the leafs of the subtree. The criteria are specified in
-      the environment as well as the subdivision method. By default surface area
-      heuristics is used.
-        
-      @param subtree root of the subtree
-
-      @return true if subdivision was performed, false if subdivision criteria
-      were already met
-  */
-  virtual ViewCellKdNode *Subdivide(const TraversalData &tdata);
-
-  /** Get the root of the tree */
-  ViewCellKdNode *GetRoot() const {
-    return mRoot;
-  }
-
-
-  AxisAlignedBox3 GetBox() const { return mBox; }
-
-  int
-  CastRay(
-          Ray &ray
-          );
-
-  const ViewCellKdTreeStatistics &GetStatistics() const {
-    return mStat;
-  }
-
-  void
-  CollectObjects(ViewCellKdNode *n, ObjectContainer &objects);
-        
-  void
-  CollectLeaves(vector<ViewCellKdLeaf *> &leaves);
-        
-  AxisAlignedBox3 GetBox(const ViewCellKdNode *node) const {
-    ViewCellKdInterior *parent = node->mParent;
-    if (parent == NULL)
-      return mBox;
-    
-    if (!node->IsLeaf())
-      return ((ViewCellKdInterior *)node)->mBox;
-    
-    AxisAlignedBox3 box(parent->mBox);
-    if (parent->mFront == node)
-      box.SetMin(parent->mAxis, parent->mPosition);
-    else
-      box.SetMax(parent->mAxis, parent->mPosition);
-    return box;
-  }
-        
-  ViewCellKdNode *
-  FindRandomNeighbor(ViewCellKdNode *n,
-                                                                                 bool onlyUnmailed
-                                                                                 );
-  
-  ViewCellKdNode *
-  ViewCellKd::GetRandomLeaf(const Plane3 &halfspace);
-
-        ViewCellKdNode *
-        GetRandomLeaf(const bool onlyUnmailed = false);
-
-  int
-  FindNeighbors(ViewCellKdNode *n,
-                vector<ViewCellKdNode *> &neighbors,
-                bool onlyUnmailed
-                );
-
-  int
-  CollectLeafPvs();
-
-protected:
-
-  struct RayData {
-    // pointer to the actual ray
-    Ray *ray;
-    
-    // endpoints  - do we need them?
-#if USE_FIXEDPOINT_T
-    short tmin, tmax;
-#else
-    float tmin, tmax;
-#endif
-
-    RayData():ray(NULL) {}
-    RayData(Ray *r):ray(r), tmin(0),
-
-#if USE_FIXEDPOINT_T
-#define FIXEDPOINT_ONE 0x7FFE
-                          //                      tmax(0xFFFF)
-                          tmax(FIXEDPOINT_ONE)
-#else
-      tmax(1.0f)
-#endif
-    {}
-
-    RayData(Ray *r,
-            const float _min,
-            const float _max
-            ):ray(r) {
-      SetTMin(_min);
-      SetTMax(_max);
-    }
-
-    RayData(Ray *r,
-            const short _min,
-            const float _max
-            ):ray(r), tmin(_min) {
-      SetTMax(_max);
-    }
-
-    RayData(Ray *r,
-            const float _min,
-            const short _max
-            ):ray(r), tmax(_max) {
-      SetTMin(_min);
-    }
-
-    friend bool operator<(const RayData &a, const RayData &b) {
-      return a.ray < b.ray;
-    }
-    
-    
-    float ExtrapOrigin(const int axis) const {
-      return ray->GetLoc(axis) + GetTMin()*ray->GetDir(axis);
-    }
-    
-    float ExtrapTermination(const int axis) const {
-      return ray->GetLoc(axis) + GetTMax()*ray->GetDir(axis);
-    }
-    
-#if USE_FIXEDPOINT_T
-    float GetTMin () const { return tmin/(float)(FIXEDPOINT_ONE); }
-    float GetTMax () const { return tmax/(float)(FIXEDPOINT_ONE); }
-
-    void SetTMin (const float t) {
-      tmin = (short) (t*(float)(FIXEDPOINT_ONE));
-    }
-    
-    void SetTMax (const float t) {
-      tmax = (short) (t*(float)(FIXEDPOINT_ONE));
-      tmax++;
-      //      if (tmax!=0xFFFF)
-      //        tmax++;
-    }
-#else
-    float GetTMin () const { return tmin; }
-    float GetTMax () const { return tmax; }
-
-    void SetTMin (const float t) { tmin = t; }
-    void SetTMax (const float t) { tmax = t; }
-#endif 
-  };
-
-  struct RayTraversalData {
-    ViewCellKdNode *mNode;
-    Vector3 mExitPoint;
-    float mMaxT;
+    VspKdTreeNode::RayInfo rayData;
+    VspKdTreeNode *node;
     
     RayTraversalData() {}
-    RayTraversalData(ViewCellKdNode *n,
-                     const Vector3 &p,
-                     const float maxt):
-      mNode(n), mExitPoint(p), mMaxT(maxt) {}
+    RayTraversalData(VspKdTreeNode *n,
+                                                                                 const VspKdTreeNode::RayInfo &data):
+      rayData(data), node(n) {}
   };
-
-  // --------------------------------------------------------------
-  // For sorting objects
-  // --------------------------------------------------------------
-  struct  SortableEntry
-  {
-    enum {
-      BOX_MIN,
-      BOX_MAX
-    };
-    
-    int type;
-    float value;
-    Intersectable *intersectable;
-    
-    SortableEntry() {}
-    SortableEntry(const int t, const float v, Intersectable *i):type(t),
-                                                                value(v),
-                                                                intersectable(i) {}
-    
-    bool operator<(const SortableEntry &b) const {
-      return value < b.value;
-    }
-    
-  };
+        
+public:
+  /////////////////////////////
+  // The core pointer
+  VspKdTreeNode *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 };
+  int splitType;
+        
+  // 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;
-
-  float
-  BestCostRatio(
-                ViewCellKdLeaf *node,
-                const AxisAlignedBox3 &box,
-                const int axis,
-                float &position,
-                int &objectsBack,
-                int &objectsFront
-                );
+  /////////////////////////////
+
+  VspKdStatistics stat;
+        
+  
+  VspKdTree();
+  virtual ~VspKdTree();
+
+  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);
+        }
+
+  
+        
+  VspKdTreeNode *
+  Locate(const Vector3 &v);
+        
+  VspKdTreeNode *
+  SubdivideNode(VspKdTreeLeaf *leaf,
+                                                                const AxisAlignedBox3 &box,
+                                                                AxisAlignedBox3 &backBox,
+                                                                AxisAlignedBox3 &frontBox
+                                                                );
+        
+  VspKdTreeNode *
+  Subdivide(const TraversalData &tdata);
+        
+  int
+  SelectPlane(VspKdTreeLeaf *leaf,
+                                                        const AxisAlignedBox3 &box,
+                                                        float &position,
+                                                        int &raysBack,
+                                                        int &raysFront,
+                                                        int &pvsBack,
+                                                        int &pvsFront
+                                                        );
 
   void
   SortSplitCandidates(
-                      ViewCellKdLeaf *node,
-                      const int axis
-                      );
+                                                                                        VspKdTreeLeaf *node,
+                                                                                        const int axis
+                                                                                        );
+        
+  
+  // return memory usage in MB
+  float GetMemUsage() const {
+    return 
+      (sizeof(VspKdTree) +
+       stat.Leaves()*sizeof(VspKdTreeLeaf) +
+       stat.Interior()*sizeof(VspKdTreeInterior) +
+       stat.rayRefs*sizeof(VspKdTreeNode::RayInfo))/(1024.0f*1024.0f);
+  }
+        
+  float GetRayMemUsage() const {
+    return 
+      stat.rays*(sizeof(VssRay))/(1024.0f*1024.0f);
+  }
+  
+        float
+        BestCostRatioHeuristic(
+                                                                                                 VspKdTreeLeaf *node,
+                                                                                                 int &axis,
+                                                                                                 float &position,
+                                                                                                 int &raysBack,
+                                                                                                 int &raysFront,
+                                                                                                 int &pvsBack,
+                                                                                                 int &pvsFront
+                                                                                                 );
+
+        float
+        BestCostRatioRegular(
+                                                                                         VspKdTreeLeaf *node,
+                                                                                         int &axis,
+                                                                                         float &position,
+                                                                                         int &raysBack,
+                                                                                         int &raysFront,
+                                                                                         int &pvsBack,
+                                                                                         int &pvsFront
+
+                                                                                         );
+        
+        float
+        EvalCostRatio(
+                                                                VspKdTreeLeaf *node,
+                                                                const int axis,
+                                                                const float position,
+                                                                int &raysBack,
+                                                                int &raysFront,
+                                                                int &pvsBack,
+                                                                int &pvsFront
+                                                                );
+
+  AxisAlignedBox3 GetBBox(const VspKdTreeNode *node) {
+    if (node->parent == NULL)
+      return bbox;
+
+    if (!node->IsLeaf())
+      return ((VspKdTreeInterior *)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 GetDirBBox(const VspKdTreeNode *node) {
+
+    if (node->parent == NULL)
+      return dirBBox;
+    
+    if (!node->IsLeaf() )
+      return ((VspKdTreeInterior *)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(VspKdTreeNode *node, const int time);
 
   void
+  CountAccess(VspKdTreeInterior *node, const long time) {
+    node->accesses++;
+    node->lastAccessTime = time;
+  }
+
+  VspKdTreeNode *
+  SubdivideLeaf(
+                                                                VspKdTreeLeaf *leaf,
+                                                                const float SAThreshold
+                                                                );
+
+  void
+  RemoveRay(VssRay *ray,
+                                                vector<VspKdTreeLeaf *> *affectedLeaves,
+                                                const bool removeAllScheduledRays
+                                                );
+
+  void
+  AddRay(VssRay *ray);
+        
+  void
+  TraverseInternalNode(
+                                                                                         RayTraversalData &data,
+                                                                                         stack<RayTraversalData> &tstack);
+
+        void
   EvaluateLeafStats(const TraversalData &data);
 
-  ViewCellKdNode *
-  SubdivideNode(
-                ViewCellKdLeaf *leaf,
-                const AxisAlignedBox3 &box,
-                AxisAlignedBox3 &backBBox,
-                AxisAlignedBox3 &frontBBox
-                );
-
-  bool
-  TerminationCriteriaMet(const ViewCellKdLeaf *leaf);
-  
-  int
-  SelectPlane(ViewCellKdLeaf *leaf,
-              const AxisAlignedBox3 &box,
-              float &position
-              );
-
-
-
-  float mSplitBorder;
-  int mTermMaxDepth;
-  int mTermMinCost;
-  float mMaxCostRatio;
-  float mCt_div_ci;
-  int mSplitMethod;
-  bool mSahUseFaces;
-  /// root of the tree
-  ViewCellKdNode *mRoot;
-  /// bounding box of the tree root
-  AxisAlignedBox3 mBox;
-  ViewCellKdTreeStatistics mStat;
-  
+
+        int
+        GetRootPvsSize() const {
+                return GetPvsSize(root, bbox);
+        }
+        
+        int
+        GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const;
+
+        void
+        GetRayContributionStatistics(
+                                                                                                                         float &minRayContribution,
+                                                                                                                         float &maxRayContribution,
+                                                                                                                         float &avgRayContribution
+                                                                                                                         );
+
+        int
+        GenerateRays(const float ratioPerLeaf,
+                                                         SimpleRayContainer &rays);
+
+        float
+        GetAvgPvsSize();
+
 };
 
-#endif
+
+#endif // __LSDS_KDTREE_H__
+