Changeset 434 for trunk/VUT/GtpVisibilityPreprocessor/src

Timestamp:

11/25/05 14:45:50 (19 years ago)

Author:

bittner

Message:

vssbsp merge, tab change

Location:

trunk/VUT/GtpVisibilityPreprocessor/src

Files:

• Environment.cpp (modified) (1 diff)
• Exporter.h (modified) (2 diffs)
• Makefile (modified) (3 diffs)
• VssPreprocessor.cpp (modified) (7 diffs)
• VssPreprocessor.h (modified) (2 diffs)
• VssRay.h (modified) (2 diffs)
• VssTree.cpp (modified) (51 diffs)
• VssTree.h (modified) (12 diffs)
• X3dExporter.cpp (modified) (3 diffs)
• X3dExporter.h (modified) (2 diffs)
• common.h (modified) (1 diff)
• default.env (modified) (4 diffs)

trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp

@@ -1449 +1449 @@
   RegisterOption("VssTree.randomize", optBool, "randomize", "false");
   RegisterOption("VssTree.splitType", optString, "split=", "queries");
+  RegisterOption("VssTree.splitUseOnlyDrivingAxis", optBool, "splitdriving=", "false");
+
+  RegisterOption("VssTree.useRss", optBool, "rss=", "false");
+
+
+
   RegisterOption("VssTree.numberOfEndPointDomains", optInt, "endpoints=", "10000");
 

trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h

@@ -19 +19 @@
 class BspLeaf;
 class Polygon3;
+class VssTree;
 
 class Exporter
@@ -47 +48 @@
   virtual bool
   ExportKdTree(const KdTree &tree) = 0;
+
+        virtual bool
+  ExportVssTree(const VssTree &tree) = 0;
 
   virtual bool

trunk/VUT/GtpVisibilityPreprocessor/src/Makefile

@@ -1 +1 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (1.07a) (Qt 3.3.2) on: Tue Nov 22 12:10:40 2005
+# Generated by qmake (1.07a) (Qt 3.3.2) on: Wed Nov 23 09:47:01 2005
 # Project:  preprocessor.pro
 # Template: app
@@ -369 +369 @@
                 VssRay.h \
                 VspKdTree.h \
+                VssTree.h \
                 Containers.h \
                 AxisAlignedBox3.h \
@@ -439 +440 @@
                 Ray.h \
                 Plane3.h \
+                VssRay.h \
                 Matrix4x4.h \
                 Vector3.h \

trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp

@@ -11 +11 @@
 
 
-bool useViewSpaceBox = false;//true;
+bool useViewSpaceBox = true;//true;
 bool use2dSampling = false;
-bool useViewspacePlane = false;//true;
+bool useViewspacePlane = true;//true;
 
 VssPreprocessor::VssPreprocessor():
-        mPass(0),
-        mVssRays()
+  mPass(0),
+  mVssRays()
 {
   // this should increase coherence of the samples
@@ -32 +32 @@
 VssPreprocessor::~VssPreprocessor()
 {
-        CLEAR_CONTAINER(mVssRays);
+  CLEAR_CONTAINER(mVssRays);
 }
 
 void
 VssPreprocessor::SetupRay(Ray &ray, 
-                                                                                                        const Vector3 &point, 
-                                                                                                        const Vector3 &direction
-                                                                                                        )
+                                                  const Vector3 &point, 
+                                                  const Vector3 &direction
+                                                  )
 {
   ray.intersections.clear();
-        // do not store anything else then intersections at the ray
+  // do not store anything else then intersections at the ray
   ray.Init(point, direction, Ray::LOCAL_RAY);
 }
@@ -48 +48 @@
 int
 VssPreprocessor::CastRay(
-                                                                                                 Vector3 &viewPoint,
-                                                                                                 Vector3 &direction,
-                                                                                                 VssRayContainer &vssRays
-                                                                                                 )
-{
-        int hits = 0;
-        static Ray ray;
-        AxisAlignedBox3 box = mKdTree->GetBox();
-
-        AxisAlignedBox3 sbox = box;
-        sbox.Enlarge(Vector3(-Limits::Small));
-        if (!sbox.IsInside(viewPoint))
-                return 0;
-        
-        SetupRay(ray, viewPoint, direction);
-        // cast ray to KD tree to find intersection with other objects
-        Intersectable *objectA, *objectB;
-        Vector3 pointA, pointB;
-        float bsize = Magnitude(box.Size());
-        if (mKdTree->CastRay(ray)) {
-                objectA = ray.intersections[0].mObject;
-                pointA = ray.Extrap(ray.intersections[0].mT);
-        } else {
-                objectA = NULL;
-                // compute intersection with the scene bounding box
-                float tmin, tmax;
-                box.ComputeMinMaxT(ray, &tmin, &tmax);
-                if (tmax > bsize) {
-                        //                      cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
-                        //                      cerr<<"ray"<<ray<<endl;
-                }
-                pointA = ray.Extrap(tmax);
-                
-        }
-
-        bool detectEmptyViewSpace = true;
-        
-        if (detectEmptyViewSpace) {
-                SetupRay(ray, pointA, -direction);
-        } else
-                SetupRay(ray, viewPoint, -direction);
-        
-        
-        if (mKdTree->CastRay(ray)) {
-                
-                objectB = ray.intersections[0].mObject;
-          pointB = ray.Extrap(ray.intersections[0].mT);
-
-        } else {
-                objectB = NULL;
-                float tmin, tmax;
-                box.ComputeMinMaxT(ray, &tmin, &tmax);
-                if (tmax > bsize) {
-                        //                      cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
-                        //                      cerr<<"ray"<<ray<<endl;
-                }
-                
-                pointB = ray.Extrap(tmax);
-        }
-
-        VssRay *vssRay  = NULL;
-
-        bool validSample = true;
-        if (detectEmptyViewSpace) {
-                if (Distance(pointA, pointB) <
-                                Distance(viewPoint, pointA) + Distance(viewPoint, pointB) - Limits::Small) {
-                        validSample = false;
-                }
-        }
-        
-        if (validSample) {
-                if (objectA) {
-                        vssRay = new VssRay(pointB,
-                                                                                                        pointA,
-                                                                                                        objectB,
-                                                                                                        objectA);
-                        vssRays.push_back(vssRay);
-                        hits ++;
-                }
-                
-                if (objectB) {
-                        vssRay = new VssRay(pointA,
-                                                                                                        pointB,
-                                                                                                        objectA,
-                                                                                                        objectB);
-                        vssRays.push_back(vssRay);
-                        hits ++;
-                }
-        }
-        
-        return hits;
+                                                 Vector3 &viewPoint,
+                                                 Vector3 &direction,
+                                                 VssRayContainer &vssRays
+                                                 )
+{
+  int hits = 0;
+  static Ray ray;
+  AxisAlignedBox3 box = mKdTree->GetBox();
+
+  AxisAlignedBox3 sbox = box;
+  sbox.Enlarge(Vector3(-Limits::Small));
+  if (!sbox.IsInside(viewPoint))
+        return 0;
+        
+  SetupRay(ray, viewPoint, direction);
+  // cast ray to KD tree to find intersection with other objects
+  Intersectable *objectA, *objectB;
+  Vector3 pointA, pointB;
+  float bsize = Magnitude(box.Size());
+  if (mKdTree->CastRay(ray)) {
+        objectA = ray.intersections[0].mObject;
+        pointA = ray.Extrap(ray.intersections[0].mT);
+  } else {
+        objectA = NULL;
+        // compute intersection with the scene bounding box
+        float tmin, tmax;
+        box.ComputeMinMaxT(ray, &tmin, &tmax);
+        if (tmax > bsize) {
+          //                    cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
+          //                    cerr<<"ray"<<ray<<endl;
+        }
+        pointA = ray.Extrap(tmax);
+                
+  }
+
+  bool detectEmptyViewSpace = true;
+        
+  if (detectEmptyViewSpace) {
+        SetupRay(ray, pointA, -direction);
+  } else
+        SetupRay(ray, viewPoint, -direction);
+        
+        
+  if (mKdTree->CastRay(ray)) {
+                
+        objectB = ray.intersections[0].mObject;
+        pointB = ray.Extrap(ray.intersections[0].mT);
+
+  } else {
+        objectB = NULL;
+        float tmin, tmax;
+        box.ComputeMinMaxT(ray, &tmin, &tmax);
+        if (tmax > bsize) {
+          //                    cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
+          //                    cerr<<"ray"<<ray<<endl;
+        }
+                
+        pointB = ray.Extrap(tmax);
+  }
+
+  VssRay *vssRay  = NULL;
+
+  bool validSample = true;
+  if (detectEmptyViewSpace) {
+        if (Distance(pointA, pointB) <
+                Distance(viewPoint, pointA) + Distance(viewPoint, pointB) - Limits::Small) {
+          validSample = false;
+        }
+  }
+        
+  if (validSample) {
+        if (objectA) {
+          vssRay = new VssRay(pointB,
+                                                  pointA,
+                                                  objectB,
+                                                  objectA);
+          vssRays.push_back(vssRay);
+          hits ++;
+        }
+                
+        if (objectB) {
+          vssRay = new VssRay(pointA,
+                                                  pointB,
+                                                  objectA,
+                                                  objectB);
+          vssRays.push_back(vssRay);
+          hits ++;
+        }
+  }
+        
+  return hits;
 }
 
@@ -145 +145 @@
 VssPreprocessor::GetViewpoint(AxisAlignedBox3 *viewSpaceBox)
 {
-        AxisAlignedBox3 box;
-        
-        if (viewSpaceBox)
-                box =*viewSpaceBox;
-        else 
-                box = mKdTree->GetBox();
-        
-        // shrink the box in the y direction
-        return box.GetRandomPoint();
+  AxisAlignedBox3 box;
+        
+  if (viewSpaceBox)
+        box =*viewSpaceBox;
+  else 
+        box = mKdTree->GetBox();
+        
+  // shrink the box in the y direction
+  return box.GetRandomPoint();
 }
 
 Vector3
 VssPreprocessor::GetDirection(const Vector3 &viewpoint,
-                                                                                                                        AxisAlignedBox3 *viewSpaceBox
-                                                                                                                        )
-{
-        Vector3 point;
-        if (!use2dSampling) {
-                Vector3 normal;
-                int i = RandomValue(0, mObjects.size()-1);
-                Intersectable *object = mObjects[i];
-                object->GetRandomSurfacePoint(point, normal);
-        } else {
-                AxisAlignedBox3 box;
-                
-                if (viewSpaceBox)
-                        box =*viewSpaceBox;
-                else 
-                        box = mKdTree->GetBox();
-                
-                point = box.GetRandomPoint();
-                point.y = viewpoint.y;
-        }
-        
-        return point - viewpoint;
+                                                          AxisAlignedBox3 *viewSpaceBox
+                                                          )
+{
+  Vector3 point;
+  if (!use2dSampling) {
+        Vector3 normal;
+        int i = RandomValue(0, mObjects.size()-1);
+        Intersectable *object = mObjects[i];
+        object->GetRandomSurfacePoint(point, normal);
+  } else {
+        AxisAlignedBox3 box;
+                
+        if (viewSpaceBox)
+          box =*viewSpaceBox;
+        else 
+          box = mKdTree->GetBox();
+                
+        point = box.GetRandomPoint();
+        point.y = viewpoint.y;
+  }
+        
+  return point - viewpoint;
 }
 
 int
 VssPreprocessor::GenerateImportanceRays(VssTree *vssTree,
-                                                                                                                                                                const int desiredSamples,
-                                                                                                                                                                SimpleRayContainer &rays
-                                                                                                                                                                )
-{
-        int num;
-        if (0) {
-                float minRayContribution;
-                float maxRayContribution;
-                float avgRayContribution;
-                
-                vssTree->GetRayContributionStatistics(minRayContribution,
-                                                                                                                                                                        maxRayContribution,
-                                                                                                                                                                        avgRayContribution);
-                
-                cout<<
-                        "#MIN_RAY_CONTRIB\n"<<minRayContribution<<endl<<
-                        "#MAX_RAY_CONTRIB\n"<<maxRayContribution<<endl<<
-                        "#AVG_RAY_CONTRIB\n"<<avgRayContribution<<endl;
-                
-                float p = desiredSamples/(float)(avgRayContribution*vssTree->stat.Leaves());
-                num = vssTree->GenerateRays(p, rays);
-        } else {
-                int leaves = vssTree->stat.Leaves()/2;
-                num = vssTree->GenerateRays(desiredSamples, leaves, rays);
-        }
-        
-        cout<<"Generated "<<num<<" rays."<<endl;
-        
-        return num;
+                                                                                const int desiredSamples,
+                                                                                SimpleRayContainer &rays
+                                                                                )
+{
+  int num;
+  if (0) {
+        float minRayContribution;
+        float maxRayContribution;
+        float avgRayContribution;
+                
+        vssTree->GetRayContributionStatistics(minRayContribution,
+                                                                                  maxRayContribution,
+                                                                                  avgRayContribution);
+                
+        cout<<
+          "#MIN_RAY_CONTRIB\n"<<minRayContribution<<endl<<
+          "#MAX_RAY_CONTRIB\n"<<maxRayContribution<<endl<<
+          "#AVG_RAY_CONTRIB\n"<<avgRayContribution<<endl;
+                
+        float p = desiredSamples/(float)(avgRayContribution*vssTree->stat.Leaves());
+        num = vssTree->GenerateRays(p, rays);
+  } else {
+        int leaves = vssTree->stat.Leaves()/2;
+        num = vssTree->GenerateRays(desiredSamples, leaves, rays);
+  }
+        
+  cout<<"Generated "<<num<<" rays."<<endl;
+        
+  return num;
 }
 
@@ -218 +218 @@
 bool
 VssPreprocessor::ExportRays(const char *filename,
-                                                                                                                const VssRayContainer &vssRays,
-                                                                                                                const int number
-                                                                                                                )
-{
-        cout<<"Exporting vss rays..."<<endl<<flush;
-        
-        float prob = number/(float)vssRays.size();
-
-
-        Exporter *exporter = NULL;
-        exporter = Exporter::GetExporter(filename);
-        //      exporter->SetWireframe();
-        //      exporter->ExportKdTree(*mKdTree);
-        exporter->SetFilled();
-        exporter->ExportScene(mSceneGraph->mRoot);
-        exporter->SetWireframe();
-
-        if (mViewSpaceBox) {
-                exporter->SetForcedMaterial(RgbColor(1,0,0));
-                exporter->ExportBox(*mViewSpaceBox);
-                exporter->ResetForcedMaterial();
-        }
-        
-        VssRayContainer rays;   for (int i=0; i < vssRays.size(); i++)
-                if (RandomValue(0,1) < prob)
-                        rays.push_back(vssRays[i]);
-
-        exporter->ExportRays(rays, RgbColor(1, 0, 0));
-        
-        delete exporter;
-
-        cout<<"done."<<endl<<flush;
-
-        return true;
-}
-
+                                                        const VssRayContainer &vssRays,
+                                                        const int number
+                                                        )
+{
+  cout<<"Exporting vss rays..."<<endl<<flush;
+        
+  float prob = number/(float)vssRays.size();
+
+
+  Exporter *exporter = NULL;
+  exporter = Exporter::GetExporter(filename);
+  //    exporter->SetWireframe();
+  //    exporter->ExportKdTree(*mKdTree);
+  exporter->SetFilled();
+  exporter->ExportScene(mSceneGraph->mRoot);
+  exporter->SetWireframe();
+
+  if (mViewSpaceBox) {
+        exporter->SetForcedMaterial(RgbColor(1,0,0));
+        exporter->ExportBox(*mViewSpaceBox);
+        exporter->ResetForcedMaterial();
+  }
+        
+  VssRayContainer rays; for (int i=0; i < vssRays.size(); i++)
+        if (RandomValue(0,1) < prob)
+          rays.push_back(vssRays[i]);
+
+  exporter->ExportRays(rays, RgbColor(1, 0, 0));
+        
+  delete exporter;
+
+  cout<<"done."<<endl<<flush;
+
+  return true;
+}
+
+
+bool
+VssPreprocessor::ExportVssTree(char *filename,
+                                                           VssTree *tree)
+{
+  Exporter *exporter = Exporter::GetExporter(filename);
+  exporter->SetFilled();
+  exporter->ExportScene(mSceneGraph->mRoot);
+  exporter->SetWireframe();
+  bool result = exporter->ExportVssTree(*tree);
+  delete exporter;
+  return result;
+}
 
 bool
 VssPreprocessor::ExportVssTreeLeaf(char *filename,
-                                                                                                                                         VssTree *tree,
-                                                                                                                                         VssTreeLeaf *leaf)
-{
-        Exporter *exporter = NULL;
-        exporter = Exporter::GetExporter(filename);
-        exporter->SetWireframe();
-        exporter->ExportKdTree(*mKdTree);
-        
-        if (mViewSpaceBox) {
-                exporter->SetForcedMaterial(RgbColor(1,0,0));
-                exporter->ExportBox(*mViewSpaceBox);
-                exporter->ResetForcedMaterial();
-        }
-        
-        exporter->SetForcedMaterial(RgbColor(0,0,1));
-        exporter->ExportBox(tree->GetBBox(leaf));
+                                                                   VssTree *tree,
+                                                                   VssTreeLeaf *leaf)
+{
+  Exporter *exporter = NULL;
+  exporter = Exporter::GetExporter(filename);
+  exporter->SetWireframe();
+  exporter->ExportKdTree(*mKdTree);
+        
+  if (mViewSpaceBox) {
+        exporter->SetForcedMaterial(RgbColor(1,0,0));
+        exporter->ExportBox(*mViewSpaceBox);
         exporter->ResetForcedMaterial();
-        
-        VssRayContainer rays[4];
-        for (int i=0; i < leaf->rays.size(); i++) {
-                int k = leaf->rays[i].GetRayClass();
-                rays[k].push_back(leaf->rays[i].mRay);
-        }
-        
-        // SOURCE RAY
-        exporter->ExportRays(rays[0], RgbColor(1, 0, 0));
-        // TERMINATION RAY
-        exporter->ExportRays(rays[1], RgbColor(1, 1, 1));
-        // PASSING_RAY
-        exporter->ExportRays(rays[2], RgbColor(1, 1, 0));
-        // CONTAINED_RAY
-        exporter->ExportRays(rays[3], RgbColor(0, 0, 1));
-
-        delete exporter;
-        return true;
+  }
+        
+  exporter->SetForcedMaterial(RgbColor(0,0,1));
+  exporter->ExportBox(tree->GetBBox(leaf));
+  exporter->ResetForcedMaterial();
+        
+  VssRayContainer rays[4];
+  for (int i=0; i < leaf->rays.size(); i++) {
+        int k = leaf->rays[i].GetRayClass();
+        rays[k].push_back(leaf->rays[i].mRay);
+  }
+        
+  // SOURCE RAY
+  exporter->ExportRays(rays[0], RgbColor(1, 0, 0));
+  // TERMINATION RAY
+  exporter->ExportRays(rays[1], RgbColor(1, 1, 1));
+  // PASSING_RAY
+  exporter->ExportRays(rays[2], RgbColor(1, 1, 0));
+  // CONTAINED_RAY
+  exporter->ExportRays(rays[3], RgbColor(0, 0, 1));
+
+  delete exporter;
+  return true;
 }
 
@@ -297 +310 @@
 VssPreprocessor::ExportVssTreeLeaves(VssTree *tree, const int number)
 {
-        vector<VssTreeLeaf *> leaves;
-        tree->CollectLeaves(leaves);
-
-        int num = 0;
-        int i;
-        float p = number / (float)leaves.size();
-        for (i=0; i < leaves.size(); i++) {
-                if (RandomValue(0,1) < p) {
-                        char filename[64];
-                        sprintf(filename, "vss-leaf-%04d.x3d", num);
-                        ExportVssTreeLeaf(filename, tree, leaves[i]);
-                        num++;
-                }
-                if (num >= number)
-                        break;
-        }
+  vector<VssTreeLeaf *> leaves;
+  tree->CollectLeaves(leaves);
+
+  int num = 0;
+  int i;
+  float p = number / (float)leaves.size();
+  for (i=0; i < leaves.size(); i++) {
+        if (RandomValue(0,1) < p) {
+          char filename[64];
+          sprintf(filename, "vss-leaf-%04d.x3d", num);
+          ExportVssTreeLeaf(filename, tree, leaves[i]);
+          num++;
+        }
+        if (num >= number)
+          break;
+  }
+}
+
+
+float
+VssPreprocessor::GetAvgPvsSize(VssTree *tree,
+                                                           const vector<AxisAlignedBox3> &viewcells
+                                                           )
+{
+  vector<AxisAlignedBox3>::const_iterator it, it_end = viewcells.end();
+
+  int sum = 0;
+  for (it = viewcells.begin(); it != it_end; ++ it)
+        sum += tree->GetPvsSize(*it);
+        
+  return sum/(float)viewcells.size();
 }
 
@@ -326 +354 @@
 
 
-        AxisAlignedBox3 *box = new AxisAlignedBox3(mKdTree->GetBox());
-
-        if (!useViewspacePlane) {
-                float size = 0.01f;
-                float s = 0.5f - size;
-                float olds = Magnitude(box->Size());
-                box->Enlarge(box->Size()*Vector3(-s));
-                Vector3 translation = Vector3(-olds*0.1f, 0, 0);
-                box->SetMin(box->Min() + translation);
-                box->SetMax(box->Max() + translation);
-        } else {
-                
-                // sample city like heights
-                box->SetMin(1, box->Min(1) + box->Size(1)*0.1);
-                box->SetMax(1, box->Min(1) + box->Size(1)*0.2);
-        }
-
-        if (use2dSampling)
-                box->SetMax(1, box->Min(1));
-        
-        if (useViewSpaceBox)
-                mViewSpaceBox = box;
-        else
-                mViewSpaceBox = NULL;
-                
-
-        VssTree *vssTree = NULL;
-
-        RayContainer bspRays;
+  AxisAlignedBox3 *box = new AxisAlignedBox3(mKdTree->GetBox());
+
+  if (!useViewspacePlane) {
+        float size = 0.01f;
+        float s = 0.5f - size;
+        float olds = Magnitude(box->Size());
+        box->Enlarge(box->Size()*Vector3(-s));
+        Vector3 translation = Vector3(-olds*0.1f, 0, 0);
+        box->SetMin(box->Min() + translation);
+        box->SetMax(box->Max() + translation);
+  } else {
+                
+        // sample city like heights
+        box->SetMin(1, box->Min(1) + box->Size(1)*0.1);
+        box->SetMax(1, box->Min(1) + box->Size(1)*0.2);
+  }
+
+  if (use2dSampling)
+        box->SetMax(1, box->Min(1));
+        
+  if (useViewSpaceBox)
+        mViewSpaceBox = box;
+  else
+        mViewSpaceBox = NULL;
+                
+
+  VssTree *vssTree = NULL;
+
+  RayContainer bspRays;
 
   while (totalSamples < mInitialSamples) {
-                int passContributingSamples = 0;
-                int passSampleContributions = 0;
-                int passSamples = 0;
-                int index = 0;
-                
-                int sampleContributions;
-                
-                int s = Min(mSamplesPerPass, mInitialSamples);
-                for (int k=0; k < s; k++) {
+        int passContributingSamples = 0;
+        int passSampleContributions = 0;
+        int passSamples = 0;
+        int index = 0;
+                
+        int sampleContributions;
+                
+        int s = Min(mSamplesPerPass, mInitialSamples);
+        for (int k=0; k < s; k++) {
                         
-                        Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
-                        Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
+          Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
+          Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
                         
-                        sampleContributions = CastRay(viewpoint, direction, mVssRays);
+          sampleContributions = CastRay(viewpoint, direction, mVssRays);
                         
                         
-                        //-- CORR matt: put block inside loop
-                        if (sampleContributions) {
-                                passContributingSamples ++;
-                                passSampleContributions += sampleContributions;
-                        }
-                        passSamples++;
-                        totalSamples++;
-                }
+          //-- CORR matt: put block inside loop
+          if (sampleContributions) {
+                passContributingSamples ++;
+                passSampleContributions += sampleContributions;
+          }
+          passSamples++;
+          totalSamples++;
+        }
     
-                mPass++;
-                
-                int pvsSize = 0;
-                float avgRayContrib = (passContributingSamples > 0) ? 
-                        passSampleContributions/(float)passContributingSamples : 0;
-                
-                cout << "#Pass " << mPass << " : t = " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
-                cout << "#TotalSamples=" << totalSamples/1000 
-                                 << "k   #SampleContributions=" << passSampleContributions << " (" 
-                                 << 100*passContributingSamples/(float)passSamples<<"%)" << " avgPVS="
-                                 << pvsSize/(float)mObjects.size() << endl 
-                                 << "avg ray contrib=" << avgRayContrib << endl;
-                
-                mStats <<
-                        "#Pass\n" <<mPass<<endl<<
-                        "#Time\n" << TimeDiff(startTime, GetTime())*1e-3 << endl<<
-                        "#TotalSamples\n" << totalSamples<< endl<<
-                        "#SampleContributions\n" << passSampleContributions << endl << 
-                        "#PContributingSamples\n"<<100*passContributingSamples/(float)passSamples<<endl <<
-                        "#AvgPVS\n"<< pvsSize/(float)mObjects.size() << endl <<
-                        "#AvgRayContrib\n" << avgRayContrib << endl;
-
-
-
-                
-        }
-        
-        cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
-        cout << "#totalRayStackSize=" << mVssRays.size() << endl <<flush;
-
-        
-        int numExportRays = 10000;
-        //int numExportRays = 0;
+        mPass++;
+                
+        int pvsSize = 0;
+        float avgRayContrib = (passContributingSamples > 0) ? 
+          passSampleContributions/(float)passContributingSamples : 0;
+                
+        cout << "#Pass " << mPass << " : t = " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+        cout << "#TotalSamples=" << totalSamples/1000 
+                 << "k   #SampleContributions=" << passSampleContributions << " (" 
+                 << 100*passContributingSamples/(float)passSamples<<"%)" << " avgPVS="
+                 << pvsSize/(float)mObjects.size() << endl 
+                 << "avg ray contrib=" << avgRayContrib << endl;
+                
+        mStats <<
+          "#Pass\n" <<mPass<<endl<<
+          "#Time\n" << TimeDiff(startTime, GetTime())*1e-3 << endl<<
+          "#TotalSamples\n" << totalSamples<< endl<<
+          "#SampleContributions\n" << passSampleContributions << endl << 
+          "#PContributingSamples\n"<<100*passContributingSamples/(float)passSamples<<endl <<
+          "#AvgPVS\n"<< pvsSize/(float)mObjects.size() << endl <<
+          "#AvgRayContrib\n" << avgRayContrib << endl;
+
+
+
+                
+  }
+        
+  cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
+  cout << "#totalRayStackSize=" << mVssRays.size() << endl <<flush;
+        
+  int numExportRays = 10000;
+  //int numExportRays = 0;
+
+  if (numExportRays) {
+        char filename[64];
+        sprintf(filename, "vss-rays-initial.x3d");
+        ExportRays(filename, mVssRays, numExportRays);
+  }
+        
+  //-- construct BSP view cells
+  if (ViewCell::sHierarchy == ViewCell::BSP) 
+        {
+          const int bspSamples = min((int)mVssRays.size(), mBspConstructionSamples);
+                
+          for (int i = 0; i < bspSamples; ++ i)
+                bspRays.push_back(new Ray(*mVssRays[i]));
+
+          //-- construct BSP tree using the samples
+          mBspTree = new BspTree(&mUnbounded);  
+
+          mBspTree->SetGenerateViewCells(true);
+          mBspTree->Construct(bspRays);
+
+          Exporter *exporter = Exporter::GetExporter("vccbsprays.x3d");
+                        
+          // export rays piercing this view cell
+          exporter->ExportRays(bspRays, 1000, RgbColor(0, 1, 0));
+
+          // cast remaining initial rays into BSP tree
+          for (int i = bspSamples; i < (int)mVssRays.size(); ++ i)
+                CastRay(*mBspTree, *mVssRays[i]);
+        }
+
+  vssTree = new VssTree;
+  // viewcells = Construct(mVssRays);
+        
+  vssTree->Construct(mVssRays, mViewSpaceBox);
+  cout<<"VssTree root PVS size = "<<vssTree->GetRootPvsSize()<<endl;
+
+  ExportVssTree("vss-tree.x3d", vssTree);
+
+  ExportVssTreeLeaves(vssTree, 10);
+
+  // viewcells->UpdatePVS(newVssRays);
+  // get viewcells as kd tree boxes
+  vector<AxisAlignedBox3> kdViewcells;
+  if (0) {
+        vector<KdLeaf *> leaves;
+        mKdTree->CollectLeaves(leaves);
+        vector<KdLeaf *>::const_iterator it;
+        int targetLeaves = 50;
+        float prob = targetLeaves/(float)leaves.size();
+        for (it = leaves.begin(); it != leaves.end(); ++it)
+          if (RandomValue(0.0f,1.0f) < prob)
+                kdViewcells.push_back(mKdTree->GetBox(*it));
+                
+        float avgPvs = GetAvgPvsSize(vssTree, kdViewcells);
+        cout<<"Initial average PVS size = "<<avgPvs<<endl;
+  }
+
+        
+  int samples = 0;
+  int pass = 0;
+  while (1) {
+        int num = mVssSamplesPerPass;
+        SimpleRayContainer rays;
+        VssRayContainer vssRays;
+                
+        if (!mUseImportanceSampling) {
+          for (int j=0; j < num; j++) {
+                Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
+                Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
+                rays.push_back(SimpleRay(viewpoint, direction));
+          }
+        } else {
+          num = GenerateImportanceRays(vssTree, num, rays);
+        }
+                
+                
+        for (int i=0; i < rays.size(); i++)
+          CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays);
+                
+        vssTree->AddRays(vssRays);
+                
+        if (0) {
+          int subdivided = vssTree->UpdateSubdivision();
+          cout<<"subdivided leafs = "<<subdivided<<endl;
+        }
+
+        float avgPvs = GetAvgPvsSize(vssTree, kdViewcells);
+        cout<<"Average PVS size = "<<avgPvs<<endl;
 
         if (numExportRays) {
-                char filename[64];
-                sprintf(filename, "vss-rays-initial.x3d");
-                ExportRays(filename, mVssRays, numExportRays);
-        }
-
-        //-- construct BSP view cells
+          char filename[64];
+          if (mUseImportanceSampling)
+                sprintf(filename, "vss-rays-i%04d.x3d", pass);
+          else
+                sprintf(filename, "vss-rays-%04d.x3d", pass);
+                        
+          ExportRays(filename, vssRays, numExportRays);
+        }
+
+        // cast rays into BSP tree
         if (ViewCell::sHierarchy == ViewCell::BSP) 
+          {
+                for (int i = 0; i < (int)vssRays.size(); ++ i)
+                  {
+                        CastRay(*mBspTree, *vssRays[i]);
+                  }
+          }
+
+        samples+=num;
+        float pvs = vssTree->GetAvgPvsSize();
+        cout<<"*****************************\n";
+        cout<<samples<<" avgPVS ="<<pvs<<endl;
+        cout<<"VssTree root PVS size = "<<vssTree->GetRootPvsSize()<<endl;
+        cout<<"*****************************\n";
+        if (samples >= mVssSamples)
+          break;
+        pass++;
+  }
+        
+  delete vssTree;
+
+  if (ViewCell::sHierarchy == ViewCell::BSP) 
         {
-                const int bspSamples = min((int)mVssRays.size(), mBspConstructionSamples);
-        
-                for (int i = 0; i < bspSamples; ++ i)
-                        bspRays.push_back(new Ray(*mVssRays[i]));
-
-                //-- construct BSP tree using the samples
-                mBspTree = new BspTree(&mUnbounded);    
-
-                mBspTree->SetGenerateViewCells(true);
-                mBspTree->Construct(bspRays);
-
-                Exporter *exporter = Exporter::GetExporter("vccbsprays.x3d");
+          Debug << mBspTree->GetStatistics();
+
+          ObjectContainer objects;
+          ExportSplits(objects, bspRays, 10000);
+          ExportBspPvs(objects, bspRays, 10000);
+
+          BspViewCellsStatistics stat;
+          mBspTree->EvaluateViewCellsStats(stat);
+          Debug << "original view cell partition:\n" << stat << endl;
+
+          // clear BSP samples
+          CLEAR_CONTAINER(bspRays);
+        }
+
+  return true;
+}
+
+void VssPreprocessor::CastRay(const BspTree &tree, const VssRay & vssRay)
+{
+  //-- cast ray to BSP tree to get intersection with view cells
+  Ray ray(vssRay);
+
+  Debug << ray << endl;
+  if (ray.intersections.empty())
+        Debug << "empty ray" << endl;
+  else 
+        Debug << "intersection: " << ray.intersections[0].mT << " " << ray.intersections[0].mObject << endl;
+
+  mBspTree->CastRay(ray);
+                                
+  //if (0 && ray.sourceObject.mObject)
+  //sampleContributions += 
+  //AddObjectSamples(ray.sourceObject.mObject, ray);
                         
-                // export rays piercing this view cell
-                exporter->ExportRays(bspRays, 1000, RgbColor(0, 1, 0));
-
-                // cast remaining initial rays into BSP tree
-                for (int i = bspSamples; i < (int)mVssRays.size(); ++ i)
-                        CastRay(*mBspTree, *mVssRays[i]);
-        }
-
-        vssTree = new VssTree;
-
-        vssTree->Construct(mVssRays, mViewSpaceBox);
-
-        cout<<"VssTree root PVS size = "<<vssTree->GetRootPvsSize()<<endl;
-
-        ExportVssTreeLeaves(vssTree, 10);
-        
-        int samples = 0;
-        int pass = 0;
-        while (1) {
-                int num = mVssSamplesPerPass;
-                SimpleRayContainer rays;
-                VssRayContainer vssRays;
-                
-                if (!mUseImportanceSampling) {
-                        for (int j=0; j < num; j++) {
-                                Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
-                                Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
-                                rays.push_back(SimpleRay(viewpoint, direction));
-                        }
-                } else {
-                        num = GenerateImportanceRays(vssTree, num, rays);
-                }
-                
-                
-                for (int i=0; i < rays.size(); i++)
-                        CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays);
-                
-                vssTree->AddRays(vssRays);
-                
-                if (1) {
-                        int subdivided = vssTree->UpdateSubdivision();
-                        cout<<"subdivided leafs = "<<subdivided<<endl;
-                }
-                
-                if (numExportRays) {
-                        char filename[64];
-                        if (mUseImportanceSampling)
-                                sprintf(filename, "vss-rays-i%04d.x3d", pass);
-                        else
-                                sprintf(filename, "vss-rays-%04d.x3d", pass);
-                        
-                        ExportRays(filename, vssRays, numExportRays);
-                }
-
-                // cast rays into BSP tree
-                if (ViewCell::sHierarchy == ViewCell::BSP) 
-                {
-                        for (int i = 0; i < (int)vssRays.size(); ++ i)
-                        {
-                                CastRay(*mBspTree, *vssRays[i]);
-                        }
-                }
-
-                samples+=num;
-                float pvs = vssTree->GetAvgPvsSize();
-                cout<<"*****************************\n";
-                cout<<samples<<" avgPVS ="<<pvs<<endl;
-                cout<<"VssTree root PVS size = "<<vssTree->GetRootPvsSize()<<endl;
-                cout<<"*****************************\n";
-                if (samples >= mVssSamples)
-                        break;
-                pass++;
-        }
-
-        delete vssTree;
-
-        if (ViewCell::sHierarchy == ViewCell::BSP) 
+  if (!ray.intersections.empty()) // second intersection found
         {
-                Debug << mBspTree->GetStatistics();
-
-                ObjectContainer objects;
-                ExportSplits(objects, bspRays, 10000);
-                ExportBspPvs(objects, bspRays, 10000);
-
-                BspViewCellsStatistics stat;
-                mBspTree->EvaluateViewCellsStats(stat);
-                Debug << "original view cell partition:\n" << stat << endl;
-
-                // clear BSP samples
-                CLEAR_CONTAINER(bspRays);
-        }
-
-        return true;
-}
-
-void VssPreprocessor::CastRay(const BspTree &tree, const VssRay & vssRay)
-{
-        //-- cast ray to BSP tree to get intersection with view cells
-        Ray ray(vssRay);
-
-        Debug << ray << endl;
-        if (ray.intersections.empty())
-                Debug << "empty ray" << endl;
-        else 
-                Debug << "intersection: " << ray.intersections[0].mT << " " << ray.intersections[0].mObject << endl;
-
-        mBspTree->CastRay(ray);
-                                
-        //if (0 && ray.sourceObject.mObject)
-                //sampleContributions += 
-                //AddObjectSamples(ray.sourceObject.mObject, ray);
-                        
-                if (!ray.intersections.empty()) // second intersection found
-                {
-                        //sampleContributions += 
-                                AddObjectSamples(ray.intersections[0].mObject, ray);
-                }
+          //sampleContributions += 
+          AddObjectSamples(ray.intersections[0].mObject, ray);
+        }
 }
 
 int VssPreprocessor::AddObjectSamples(Intersectable *obj, const Ray &ray)
 {
-        int contributingSamples = 0;
-        int j;
+  int contributingSamples = 0;
+  int j;
  
-        // object can be seen from the view cell => add to view cell pvs
-        for (j=0; j < ray.bspIntersections.size(); ++ j)
+  // object can be seen from the view cell => add to view cell pvs
+  for (j=0; j < ray.bspIntersections.size(); ++ j)
         {       
+          BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
+          // if ray not in unbounded space
+          if (leaf->GetViewCell() != &mUnbounded)
+                contributingSamples += 
+                  leaf->GetViewCell()->GetPvs().AddSample(obj);
+        }
+ 
+  // rays passing through this viewcell
+  if (mPass > 1)
+        for (j=1; j < ((int)ray.bspIntersections.size() - 1); ++ j) 
+          {
                 BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
-                // if ray not in unbounded space
+                
                 if (leaf->GetViewCell() != &mUnbounded)
-                        contributingSamples += 
-                                leaf->GetViewCell()->GetPvs().AddSample(obj);
-        }
- 
-        // rays passing through this viewcell
-        if (mPass > 1)
-                for (j=1; j < ((int)ray.bspIntersections.size() - 1); ++ j) 
-                {
-                        BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
-
-                        if (leaf->GetViewCell() != &mUnbounded)
-                                leaf->GetViewCell()->
-                                        AddPassingRay(ray, contributingSamples ? 1 : 0);
-                }
- 
-        return contributingSamples;
-}
+                  leaf->GetViewCell()->
+                        AddPassingRay(ray, contributingSamples ? 1 : 0);
+          }
+  
+  return contributingSamples;
+}

trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.h

@@ -14 +14 @@
 class VssPreprocessor : public Preprocessor {
 public:
-        int mPass;
+  int mPass;
   int mSamplesPerPass;
   int mVssSamplesPerPass;
   int mInitialSamples;
   int mVssSamples;
-        bool mUseImportanceSampling;
+  bool mUseImportanceSampling;
 
-        AxisAlignedBox3 *mViewSpaceBox;
+  AxisAlignedBox3 *mViewSpaceBox;
 
-        ofstream mStats;
+  ofstream mStats;
         
   ObjectContainer mObjects;
 
-        // rays cast during the processing
+  // rays cast during the processing
   VssRayContainer mVssRays;
         
@@ -35 +35 @@
   virtual bool ComputeVisibility();
 
-        Vector3
-        GetViewpoint(AxisAlignedBox3 *viewSpaceBox);
+  Vector3
+  GetViewpoint(AxisAlignedBox3 *viewSpaceBox);
 
-        Vector3
-        GetDirection(const Vector3 &viewpoint,
-                                                         AxisAlignedBox3 *viewSpaceBox
-                                                         );
+  Vector3
+  GetDirection(const Vector3 &viewpoint,
+                           AxisAlignedBox3 *viewSpaceBox
+                           );
         
   void
   SetupRay(Ray &ray, 
-                                         const Vector3 &point, 
-                                         const Vector3 &direction
-                                         );
+                   const Vector3 &point, 
+                   const Vector3 &direction
+                   );
 
 
         
-        int
-        CastRay(
-                                        Vector3 &viewPoint,
-                                        Vector3 &direction,
-                                        VssRayContainer &vssRays
+  int
+  CastRay(
+                  Vector3 &viewPoint,
+                  Vector3 &direction,
+                  VssRayContainer &vssRays
 
-                                        );
+                  );
                 
 
         
-        virtual bool BuildBspTree() { return false; }
+  virtual bool BuildBspTree() { return false; }
 
 
-        bool
-        ExportRays(const char *filename,
-                                                 const VssRayContainer &vssRays,
-                                                 const int number
+  bool
+  ExportRays(const char *filename,
+                         const VssRayContainer &vssRays,
+                         const int number
+                         );
+
+  int
+  GenerateImportanceRays(VssTree *vssTree,
+                                                 const int desiredSamples,
+                                                 SimpleRayContainer &rays
                                                  );
 
-        int
-        GenerateImportanceRays(VssTree *vssTree,
-                                                                                                 const int desiredSamples,
-                                                                                                 SimpleRayContainer &rays
-                                                                                                 );
+
+  bool
+  ExportVssTreeLeaf(char *filename,
+                                        VssTree *tree,
+                                        VssTreeLeaf *leaf);
+
+  void
+  ExportVssTreeLeaves(VssTree *tree, const int number);
+        
+
+  bool
+  ExportVssTree(char *filename,
+                                VssTree *tree);
 
 
-        bool
-        ExportVssTreeLeaf(char *filename,
-                                                                                VssTree *tree,
-                                                                                VssTreeLeaf *leaf);
-
-        void
-        ExportVssTreeLeaves(VssTree *tree, const int number);
+  float
+  GetAvgPvsSize(VssTree *tree,
+                                const vector<AxisAlignedBox3> &viewcells
+                                );
+  
+  void CastRay(const BspTree &tree, const VssRay & vssRay);
         
-        void CastRay(const BspTree &tree, const VssRay & vssRay);
-        
-        int AddObjectSamples(Intersectable *obj, const Ray &ray);
+  int AddObjectSamples(Intersectable *obj, const Ray &ray);
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h

@@ -136 +136 @@
   Vector3 GetNormalizedDir() const { return (mTermination - mOrigin)*mInvSize; }
 
+        Vector3 Extrap(const float t) const {
+                return GetOrigin() + t * GetDir();
+        }
+        
         float GetDirParametrization(const int axis) const;
         
@@ -200 +204 @@
       }
   }
+
+        static Vector3
+        GetDirection(const float a, const float b) {
+                return Vector3(sin(a), sin(b), cos(a));
+        }
+
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp

@@ -26 +26 @@
 
 
-#define DEBUG_DIR_SPLIT 0
-
-
-
 // Static variables
 int
@@ -36 +32 @@
 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);
-                }
-        }
-}
-
+                          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
@@ -70 +65 @@
 {
   environment->GetIntValue("VssTree.maxDepth", termMaxDepth);
-        environment->GetIntValue("VssTree.minPvs", termMinPvs);
-        environment->GetIntValue("VssTree.minRays", termMinRays);
-        environment->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution);
+  environment->GetIntValue("VssTree.minPvs", termMinPvs);
+  environment->GetIntValue("VssTree.minRays", termMinRays);
+  environment->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution);
   environment->GetFloatValue("VssTree.maxCostRatio", termMaxCostRatio);
 
@@ -92 +87 @@
   float refDirAngle;
   environment->GetFloatValue("VssTree.refDirAngle", refDirAngle);
-
+  
   environment->GetIntValue("VssTree.accessTimeThreshold", accessTimeThreshold);
   //= 1000;
@@ -99 +94 @@
 
   //  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);
+  //  cosRefDir = cos(M_PI*refDirAngle/180.0);
+  //  sinRefDir = sin(M_PI*refDirAngle/180.0);
   
   
   // split type
-        char sname[128];
-        environment->GetStringValue("VssTree.splitType", sname);
+  char sname[128];
+  environment->GetStringValue("VssTree.splitType", sname);
   string name(sname);
         
@@ -113 +108 @@
     if (name.compare("heuristic") == 0)
       splitType = ESplitHeuristic;
-                else {
-                        cerr<<"Invalid VssTree split type "<<name<<endl;
-                        exit(1);
-                }
+        else
+          if (name.compare("hybrid") == 0)
+                splitType = ESplitHybrid;
+          else {
+                cerr<<"Invalid VssTree split type "<<name<<endl;
+                exit(1);
+          }
         
   environment->GetBoolValue("VssTree.randomize", randomize);
-
+  environment->GetBoolValue("VssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
+  environment->GetBoolValue("VssTree.useRss", mUseRss);
+        
   root = NULL;
   
@@ -143 +143 @@
       << rays <<endl;
 
-        app << "#N_INITPVS ( Initial PVS size )\n"
+  app << "#N_INITPVS ( Initial PVS size )\n"
       << initialPvsSize <<endl;
   
@@ -179 +179 @@
     minRaysNodes*100/(double)Leaves()<<endl;
         
-        app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
+  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"<<
+  app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
     maxRayContribNodes*100/(double)Leaves()<<endl;
 
-        app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
+  app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
     maxCostRatioNodes*100/(double)Leaves()<<endl;
 
@@ -207 +207 @@
 VssTreeLeaf::UpdatePvsSize()
 {
-        if (!mValidPvs) {
-                Intersectable::NewMail();
-                int pvsSize = 0;
-                for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
-                                ri != rays.end();
-                                ri++)
-                        if ((*ri).mRay->IsActive()) {
-                                Intersectable *object;
+  if (!mValidPvs) {
+        Intersectable::NewMail();
+        int pvsSize = 0;
+        for(VssTreeNode::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();
-                                }
+                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;
-        }
+                object = (*ri).mRay->mTerminationObject;
+                if (object && !object->Mailed()) {
+                  pvsSize++;
+                  object->Mail();
+                }
+          }
+        mPvsSize = pvsSize;
+        mValidPvs = true;
+  }
 }
 
 bool
 VssTree::ClipRay(
-                                                                 VssTreeNode::RayInfo &rayInfo,
-                                                                 const AxisAlignedBox3 &box
-                                                                 )
-{
-        float tmin, tmax;
-        static Ray ray;
-        ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
-        box.ComputeMinMaxT(ray, &tmin, &tmax);
-        if (tmin >= tmax)
-                return false;
-        
-        if (tmin > 0.0f)
-                rayInfo.SetMinT(tmin);
-
-        if (tmax < 1.0f)
-                rayInfo.SetMaxT(tmax);
-        
-        //      vssRay->SetupEndPoints(
-        //                                                                                               origin,
-        //                                                                                               termination
-        //                                                                                               );
-        return true;
+                                 VssTreeNode::RayInfo &rayInfo,
+                                 const AxisAlignedBox3 &box
+                                 )
+{
+  float tmin, tmax;
+  static Ray ray;
+  ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
+  box.ComputeMinMaxT(ray, &tmin, &tmax);
+  if (tmin >= tmax)
+        return false;
+        
+  if (tmin > 0.0f)
+        rayInfo.SetMinT(tmin);
+
+  if (tmax < 1.0f)
+        rayInfo.SetMaxT(tmax);
+        
+  //    vssRay->SetupEndPoints(
+  //                                                                                             origin,
+  //                                                                                             termination
+  //                                                                                             );
+  return true;
 }
 
@@ -262 +262 @@
 void
 VssTree::Construct(
-                                                                         VssRayContainer &rays,
-                                                                         AxisAlignedBox3 *forcedBoundingBox
-                                                                         )
+                                   VssRayContainer &rays,
+                                   AxisAlignedBox3 *forcedBoundingBox
+                                   )
 {
   stat.Start();
   
-        maxMemory = maxStaticMemory;
+  maxMemory = maxStaticMemory;
 
   if (root)
     delete root;
-
+        
   root = new VssTreeLeaf(NULL, rays.size());
   // first construct a leaf that will get subdivide
@@ -281 +281 @@
   bbox.Initialize();
   dirBBox.Initialize();
-  
+
+  if (mUseRss)
+        forcedBoundingBox = NULL;
+        
   for(VssRayContainer::const_iterator ri = rays.begin();
       ri != rays.end();
       ri++) {
-
-                VssTreeNode::RayInfo info(*ri);
-                if (forcedBoundingBox)
-                        if (!ClipRay(info, *forcedBoundingBox))
-                                continue;
-                leaf->AddRay(info);
+                
+        VssTreeNode::RayInfo info(*ri);
+        if (forcedBoundingBox)
+          if (!ClipRay(info, *forcedBoundingBox))
+                continue;
+        leaf->AddRay(info);
                 
     bbox.Include((*ri)->GetOrigin());
@@ -296 +299 @@
     
                 
-                dirBBox.Include(Vector3(
-                                                                                                                (*ri)->GetDirParametrization(0),
-                                                                                                                (*ri)->GetDirParametrization(1),
-                                                                                                                0
-                                                                                                                )
-                                                                                );
-        }
-        
-        
-        if ( forcedBoundingBox ) 
-                bbox = *forcedBoundingBox;
-        
-        cout<<"Bbox = "<<bbox<<endl;
-        cout<<"Dirr Bbox = "<<dirBBox<<endl;
+        dirBBox.Include(Vector3(
+                                                        (*ri)->GetDirParametrization(0),
+                                                        (*ri)->GetDirParametrization(1),
+                                                        0
+                                                        )
+                                        );
+  }
+        
+        
+  if ( forcedBoundingBox ) 
+        bbox = *forcedBoundingBox;
+        
+  cout<<"Bbox = "<<bbox<<endl;
+  cout<<"Dirr Bbox = "<<dirBBox<<endl;
 
   stat.rays = leaf->rays.size();
-        leaf->UpdatePvsSize();
+  leaf->UpdatePvsSize();
   stat.initialPvsSize = leaf->GetPvsSize();
   // Subdivide();
@@ -325 +328 @@
   stat.Stop();
 
-        stat.Print(cout);
-        cout<<"#Total memory="<<GetMemUsage()<<endl;
+  stat.Print(cout);
+  cout<<"#Total memory="<<GetMemUsage()<<endl;
 
 }
@@ -333 +336 @@
 VssTree::UpdateSubdivision()
 {
-        priority_queue<TraversalData> tStack;
+  priority_queue<TraversalData> tStack;
   //  stack<TraversalData> tStack;
   
-        tStack.push(TraversalData(root, bbox, 0));
+  tStack.push(TraversalData(root, bbox, 0));
         
   AxisAlignedBox3 backBox;
   AxisAlignedBox3 frontBox;
 
-        maxMemory = maxTotalMemory;
-        int subdivided = 0;
-        int lastMem = 0;
+  maxMemory = maxTotalMemory;
+  int subdivided = 0;
+  int lastMem = 0;
   while (!tStack.empty()) {
                 
-                float mem = GetMemUsage();
-                
-                if ( lastMem/10 != ((int)mem)/10) {
-                        cout<<mem<<" MB"<<endl;
-                }
-                lastMem = (int)mem;
-                
-                if (  mem > maxMemory ) {
+        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;
@@ -365 +368 @@
     tStack.pop();
 
-                if (data.node->IsLeaf()) {
-                        VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
-                                                                                                                                                                data.bbox,
-                                                                                                                                                                backBox,
-                                                                                                                                                                frontBox
-                                                                                                                                                                );
-                        if (!node->IsLeaf()) {
-                                subdivided++;
-                                VssTreeInterior *interior = (VssTreeInterior *) node;
-                                // push the children on the stack
-                                tStack.push(TraversalData(interior->back, backBox, data.depth+1));
-                                tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
-                        } else {
-                                //      EvaluateLeafStats(data);
-                        }
-                } else {
-                        VssTreeInterior *interior = (VssTreeInterior *) data.node;
-                        tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
-                        tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
-                }
-  }
-        return subdivided;
+        if (data.node->IsLeaf()) {
+          VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
+                                                                                data.bbox,
+                                                                                backBox,
+                                                                                frontBox
+                                                                                );
+          if (!node->IsLeaf()) {
+                subdivided++;
+                VssTreeInterior *interior = (VssTreeInterior *) node;
+                // push the children on the stack
+                tStack.push(TraversalData(interior->back, backBox, data.depth+1));
+                tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
+          } else {
+                //      EvaluateLeafStats(data);
+          }
+        } else {
+          VssTreeInterior *interior = (VssTreeInterior *) data.node;
+          tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
+          tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
+        }
+  }
+  return subdivided;
 }
 
@@ -404 +407 @@
 
   
-        int lastMem = 0;
+  int lastMem = 0;
   while (!tStack.empty()) {
 
-                float mem = GetMemUsage();
-                
-                if ( lastMem/10 != ((int)mem)/10) {
-                        cout<<mem<<" MB"<<endl;
-                }
-                lastMem = (int)mem;
-                
-                if (  mem > maxMemory ) {
+        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;
@@ -427 +430 @@
     
     VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
-                                                                                                                                                        data.bbox,
-                                                                                                                                                        backBox,
-                                                                                                                                                        frontBox
-                                                                                                                                                        );
+                                                                          data.bbox,
+                                                                          backBox,
+                                                                          frontBox
+                                                                          );
     if (result == NULL)
       result = node;
@@ -453 +456 @@
 int
 VssTree::SelectPlane(
-                                                                                 VssTreeLeaf *leaf,
-                                                                                 const AxisAlignedBox3 &box,
-                                                                                 float &position,
-                                                                                 int &raysBack,
-                                                                                 int &raysFront,
-                                                                                 int &pvsBack,
-                                                                                 int &pvsFront
-                                                                                 )
+                                         VssTreeLeaf *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 {
+  float costRatio;
+        
+  costRatio = BestCostRatio(leaf,
+                                                        axis,
+                                                        position,
+                                                        raysBack,
+                                                        raysFront,
+                                                        pvsBack,
+                                                        pvsFront
+                                                        );
+
+  //    cout<<axis<<" r="<<costRatio<<endl;
+        
+  if (costRatio > termMaxCostRatio) {
+        //              cout<<"Too big cost ratio "<<costRatio<<endl;
+        stat.maxCostRatioNodes++;
+        return -1;
+  }
+        
+#if 0   
+  cout<<
+        "pvs="<<leaf->mPvsSize<<
+        " rays="<<leaf->rays.size()<<
+        " rc="<<leaf->GetAvgRayContribution()<<
+        " axis="<<axis<<endl;
+#endif
+        
+  return axis;
+}
+
+
+float
+VssTree::GetCostRatio(
+                                          VssTreeLeaf *leaf,
+                                          const int axis,
+                                          const float position,
+                                          const int raysBack,
+                                          const int raysFront,
+                                          const int pvsBack,
+                                          const int pvsFront
+                                          )
+{
+  bool costImportance = true;
+
+  float ratio;
+  AxisAlignedBox3 box;
+  float minBox, maxBox;
+
+  if (axis < 3) {
+        box = GetBBox(leaf);
+        minBox = box.Min(axis);
+        maxBox = box.Max(axis);
+  }     else {
+        box = GetDirBBox(leaf);
+        minBox = box.Min(axis-3);
+        maxBox = box.Max(axis-3);
+  }
+        
+  float sizeBox = maxBox - minBox;
+
+  int pvsSize = leaf->GetPvsSize();
+
+  if (!costImportance) {
+        //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+        float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
+        float newCost = ct_div_ci + sum/sizeBox;
+        float oldCost = pvsSize;
+        ratio = newCost/oldCost;
+  } else {
+        // importance based cost
+#if 0
+        float newContrib =
+          ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
+           (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
+                
+        //                      float newContrib =
+        //                              sqr(pvsBack/(raysBack + Limits::Small)) +
+        //                              sqr(pvsFront/(raysFront + Limits::Small));
+        float oldContrib = sqr(leaf->GetAvgRayContribution());
+        ratio = oldContrib/newContrib;
+#else
+#if 1
+        float newCost = raysBack*pvsBack  + raysFront*pvsFront;
+        float oldCost = leaf->rays.size()*pvsSize;
+        ratio = newCost/oldCost;
+#else
+        float newCost = (pvsBack  + pvsFront)*0.5f;
+        float oldCost = pvsSize;
+        ratio = newCost/oldCost;
+#endif
+#endif
+  }
+        
+  return ratio;
+}
+                                                        
+
+float
+VssTree::EvalCostRatio(
+                                           VssTreeLeaf *leaf,
+                                           const int axis,
+                                           const float position,
+                                           int &raysBack,
+                                           int &raysFront,
+                                           int &pvsBack,
+                                           int &pvsFront
+                                           )
+{
+  raysBack = 0;
+  raysFront = 0;
+  pvsFront = 0;
+  pvsBack = 0;
+
+        
+  Intersectable::NewMail(3);
+        
+  if (axis <= VssTreeNode::SPLIT_Z) {
+    // this is the main ray classification loop!
+    for(VssTreeNode::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);
+      }
+                
+  } else {
+                
+        // directional split
+    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                ri != leaf->rays.end();
+                ri++)
+      if ((*ri).mRay->IsActive()) {
+                                
+                // determine the side of this ray with respect to the plane
+                int side;
+                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
+                  side = 1;
+                else
+                  side = -1;
+
+                if (side <= 0)
+                  raysBack++;
+                                
+                if (side >= 0)
+                  raysFront++;
+
+                //                              (*ri).mRay->mSide = side;
+                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+                                
+      }
+  }
+
+  float ratio = GetCostRatio(
+                                                         leaf,
+                                                         axis,
+                                                         position,
+                                                         raysBack,
+                                                         raysFront,
+                                                         pvsBack,
+                                                         pvsFront);
+
+  //    cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
+  //  float oldCost = leaf->rays.size();
+
+  //    cout<<"ratio="<<ratio<<endl;
+
+  return ratio;
+}
+
+float
+VssTree::BestCostRatio(
+                                           VssTreeLeaf *leaf,
+                                           int &axis,
+                                           float &position,
+                                           int &raysBack,
+                                           int &raysFront,
+                                           int &pvsBack,
+                                           int &pvsFront
+                                           )
+{
+  int nRaysBack[6], nRaysFront[6];
+  int nPvsBack[6], nPvsFront[6];
+  float nPosition[6];
+  float nCostRatio[6];
+  int bestAxis = -1;
+        
+  AxisAlignedBox3 sBox = GetBBox(leaf);
+  AxisAlignedBox3 dBox = GetDirBBox(leaf);
+  // int sAxis = box.Size().DrivingAxis();
+  int sAxis = sBox.Size().DrivingAxis();
+  int dAxis = dBox.Size().DrivingAxis() + 3;
+
+
+  float dirSplitBoxSize = 0.01f;
+  bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
+                
+
+  for (axis = 0; axis < 5; axis++) {
+        if (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis) {
+
+
+          if (splitType == ESplitRegular) {
+                if (axis < 3)
+                  nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
+                else
+                  nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
+                                
+                nCostRatio[axis] = EvalCostRatio(leaf,
+                                                                                 axis,
+                                                                                 nPosition[axis],
+                                                                                 nRaysBack[axis],
+                                                                                 nRaysFront[axis],
+                                                                                 nPvsBack[axis],
+                                                                                 nPvsFront[axis]
+                                                                                 );
+          } else
+                if (splitType == ESplitHeuristic) {
+                  nCostRatio[axis] = EvalCostRatioHeuristic(
+                                                                                                        leaf,
+                                                                                                        axis,
+                                                                                                        nPosition[axis],
+                                                                                                        nRaysBack[axis],
+                                                                                                        nRaysFront[axis],
+                                                                                                        nPvsBack[axis],
+                                                                                                        nPvsFront[axis]);
+                } else
+                  if (splitType == ESplitHybrid) {
+                        if (leaf->depth > 7)
+                          nCostRatio[axis] = EvalCostRatioHeuristic(
+                                                                                                                leaf,
+                                                                                                                axis,
+                                                                                                                nPosition[axis],
+                                                                                                                nRaysBack[axis],
+                                                                                                                nRaysFront[axis],
+                                                                                                                nPvsBack[axis],
+                                                                                                                nPvsFront[axis]);
+                        else {
+                          if (axis < 3)
+                                nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
+                          else
+                                nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
+                                                        
+                          nCostRatio[axis] = EvalCostRatio(leaf,
+                                                                                           axis,
+                                                                                           nPosition[axis],
+                                                                                           nRaysBack[axis],
+                                                                                           nRaysFront[axis],
+                                                                                           nPvsBack[axis],
+                                                                                           nPvsFront[axis]
+                                                                                           );
+                        }
+                  } else {
                         cerr<<"VssTree: Unknown split heuristics\n";
                         exit(1);
-                }
-  }
-
-        if (costRatio > termMaxCostRatio) {
-                //              cout<<"Too big cost ratio "<<costRatio<<endl;
-                stat.maxCostRatioNodes++;
-                return -1;
-        }
-
-#if 0   
-        cout<<
-                "pvs="<<leaf->mPvsSize<<
-                " rays="<<leaf->rays.size()<<
-                " rc="<<leaf->GetAvgRayContribution()<<
-                " axis="<<axis<<endl;
-#endif
-        
-        return axis;
-}
-
-
-                                                        
-
+                  }
+                        
+                                
+          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
-VssTree::EvalCostRatio(
-                                                                                         VssTreeLeaf *leaf,
-                                                                                         const int axis,
-                                                                                         const float position,
-                                                                                         int &raysBack,
-                                                                                         int &raysFront,
-                                                                                         int &pvsBack,
-                                                                                         int &pvsFront
-                                                                                         )
-{
-        raysBack = 0;
-        raysFront = 0;
-        pvsFront = 0;
-        pvsBack = 0;
-
-
-        Intersectable::NewMail(3);
-        bool costImportance = true;
-        // eval pvs size
-        int pvsSize = leaf->GetPvsSize();
-        float ratio;
-
-        Intersectable::NewMail(3);
-
-        if (axis <= VssTreeNode::SPLIT_Z) {
-    // this is the main ray classification loop!
-    for(VssTreeNode::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;
-                
-                if (!costImportance) {
+VssTree::EvalCostRatioHeuristic(
+                                                                VssTreeLeaf *leaf,
+                                                                const int axis,
+                                                                float &bestPosition,
+                                                                int &raysBack,
+                                                                int &raysFront,
+                                                                int &pvsBack,
+                                                                int &pvsFront
+                                                                )
+{
+  AxisAlignedBox3 box;
+  float minBox, maxBox;
+        
+  if (axis < 3) {
+        box = GetBBox(leaf);
+        minBox = box.Min(axis);
+        maxBox = box.Max(axis);
+  }
+  else {
+        box = GetDirBBox(leaf);
+        minBox = box.Min(axis-3);
+        maxBox = box.Max(axis-3);
+  }
+        
+  SortSplitCandidates(leaf, axis);
+  
+  // go through the lists, count the number of objects left and right
+  // and evaluate the following cost funcion:
+  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
+        
+  int rl=0, rr = leaf->rays.size();
+  int pl=0, pr = leaf->GetPvsSize();
+  float sizeBox = maxBox - minBox;
+        
+  float minBand = minBox + 0.1*(maxBox - minBox);
+  float maxBand = minBox + 0.9*(maxBox - minBox);
+        
+  float minRatio = 1e20;
+        
+  Intersectable::NewMail();
+  // set all object as belonging to the fron pvs
+  for(VssTreeNode::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;
+        }
+        }
+
+        float position = (*ci).value;
+                
+        if (position > minBand && position < maxBand) {
                         
-                        //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
-                        float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
-                        float newCost = ct_div_ci + sum/sizeBox;
-                        float oldCost = pvsSize;
-                        ratio = newCost/oldCost;
-                } else {
-                        // importance based cost
-#if 0
-                        float newContrib =
-                                ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
-                                 (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
-
-                        //                      float newContrib =
-                        //                              sqr(pvsBack/(raysBack + Limits::Small)) +
-                        //                              sqr(pvsFront/(raysFront + Limits::Small));
-                        float oldContrib = sqr(leaf->GetAvgRayContribution());
-                        ratio = oldContrib/newContrib;
-#else
-                        float newCost = raysBack*pvsBack  + raysFront*pvsFront;
-                        float oldCost = leaf->rays.size()*pvsSize;
-                        ratio = newCost/oldCost;
-#endif
-                }
-                
-  } else {
-                
-                // directional split
-    for(VssTreeNode::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;
-
-                if (!costImportance) {
+          float ratio = GetCostRatio(
+                                                                 leaf,
+                                                                 axis,
+                                                                 position,
+                                                                 rl,
+                                                                 rr,
+                                                                 pl,
+                                                                 pr);
                         
-                        float directionalPenalty = 1.0f;
-                        float sum =
-                                (pvsBack*(position - minBox) +
-                                 pvsFront*(maxBox - position))*directionalPenalty;
-                        float newCost = ct_div_ci + sum/sizeBox;
-                        float oldCost = pvsSize;
-                        ratio = newCost/oldCost;
-                } else {
-#if 0
-                        float directionalPenalty = 0.7f;
-                        // importance based cost
-                        float newContrib =
-                                directionalPenalty*((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
-                                                                                                                (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
-                        float oldContrib = sqr(leaf->GetAvgRayContribution());
-                        ratio = oldContrib/newContrib;
-#else
-                        float newCost = raysBack*pvsBack  + raysFront*pvsFront;
-                        float oldCost = leaf->rays.size()*pvsSize;
-                        ratio = newCost/oldCost;
-#endif
-                }
-        }
-        //      cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
-        //  float oldCost = leaf->rays.size();
-
-        //      cout<<"ratio="<<ratio<<endl;
-        return ratio;
-}
-
-float
-VssTree::BestCostRatioRegular(
-                                                                                                                        VssTreeLeaf *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 = false; 
-
-        float dirSplitBoxSize = 0.01f;
-        bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
-                
-
-        for (axis = 0; axis < 5; 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]
-                                                                                                                                                         );
+          //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
+          //      cout<<"cost= "<<sum<<endl;
                         
-                        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
-VssTree::BestCostRatioHeuristic(
-                                                                                                                                VssTreeLeaf *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(VssTreeNode::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;
+          if (ratio < minRatio) {
+                minRatio = ratio;
+                bestPosition = position;
+                                
+                raysBack = rl;
+                raysFront = rr;
+                                
+                pvsBack = pl;
+                pvsFront = pr;
                                 
       }
     }
   }
-  
-  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;
+  return minRatio;
 }
 
 void
 VssTree::SortSplitCandidates(
-                                                                                                                 VssTreeLeaf *node,
-                                                                                                                 const int axis
-                                                                                                                 )
+                                                         VssTreeLeaf *node,
+                                                         const int axis
+                                                         )
 {
   
@@ -851 +903 @@
       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)
-                                                                                                                         );
-  }
-        
+        if ((*ri).mRay->IsActive()) {
+          if (axis < 3) {
+                bool positive = (*ri).mRay->HasPosDir(axis);
+                splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
+                                                                                                 SortableEntry::ERayMax,
+                                                                                                 (*ri).ExtrapOrigin(axis),
+                                                                                                 (void *)(*ri).mRay)
+                                                                   );
+                                
+                splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
+                                                                                                 SortableEntry::ERayMin,
+                                                                                                 (*ri).ExtrapTermination(axis),
+                                                                                                 (void *)(*ri).mRay)
+                                                                   );
+          } else {
+                float pos = (*ri).mRay->GetDirParametrization(axis-3);
+                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
+                                                                                                 pos - Limits::Small,
+                                                                                                 (void *)(*ri).mRay)
+                                                                   );
+                                
+                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMax,
+                                                                                                 pos + Limits::Small,
+                                                                                                 (void *)(*ri).mRay)
+                                                                   );
+          }
+        }
+  }
+
   stable_sort(splitCandidates->begin(), splitCandidates->end());
 }
@@ -879 +946 @@
     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)
+  //  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();
 
@@ -902 +969 @@
 VssTree::TerminationCriteriaSatisfied(VssTreeLeaf *leaf)
 {
-        return ( (leaf->GetPvsSize() < termMinPvs) ||
-                                        (leaf->rays.size() < termMinRays) ||
-                                        //                       (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
-                                        (leaf->depth >= termMaxDepth) ||
-                                        SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize );
+  return ( (leaf->GetPvsSize() < termMinPvs) ||
+                   (leaf->rays.size() < termMinRays) ||
+                   //                    (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
+                   (leaf->depth >= termMaxDepth) ||
+                   (SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize) ||
+                   (mUseRss && leaf->mPassingRays == leaf->rays.size())
+                   );
 }
 
@@ -912 +981 @@
 VssTreeNode *
 VssTree::SubdivideNode(
-                                                                                         VssTreeLeaf *leaf,
-                                                                                         const AxisAlignedBox3 &box,
-                                                                                         AxisAlignedBox3 &backBBox,
-                                                                                         AxisAlignedBox3 &frontBBox
-                                                                                         )
-{
-  
-        if (TerminationCriteriaSatisfied(leaf)) {
+                                           VssTreeLeaf *leaf,
+                                           const AxisAlignedBox3 &box,
+                                           AxisAlignedBox3 &backBBox,
+                                           AxisAlignedBox3 &frontBBox
+                                           )
+{
+  
+  if (TerminationCriteriaSatisfied(leaf)) {
 #if 0
-                if (leaf->depth >= termMaxDepth) {
-                        cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
-                        cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
-                }
+        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
+        return leaf;
+  }
+        
+  float position;
+        
+  // first count ray sides
   int raysBack;
   int raysFront;
-        int pvsBack;
-        int pvsFront;
+  int pvsBack;
+  int pvsFront;
         
   // select subdivision axis
   int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
-        //      cout<<axis<<" ";
-        
-        //      cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
-        //              pvsFront<<endl;
+  //    cout<<axis<<" ";
+        
+  //    cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
+  //            pvsFront<<endl;
 
   if (axis == -1) {
@@ -964 +1033 @@
 
   VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
-        back->SetPvsSize(pvsBack);
   VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
-        front->SetPvsSize(pvsFront);
 
   // replace a link from node's parent
@@ -975 +1042 @@
         
   if (axis <= VssTreeNode::SPLIT_Z) {
-                backBBox.SetMax(axis, position);
+        backBBox.SetMax(axis, position);
     frontBBox.SetMin(axis, position);
                 
-                for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                                ri != leaf->rays.end();
-                                ri++) {
+        for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                ri != leaf->rays.end();
+                ri++) {
       if ((*ri).mRay->IsActive()) {
                                 
-                                // first unref ray from the former leaf
-                                (*ri).mRay->Unref();
-
-                                Debug << "computed t: " << (*ri).mRay->mT << endl;
-                                // determine the side of this ray with respect to the plane
-                                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
-
-                                if (side == 0) {
-                                        if ((*ri).mRay->HasPosDir(axis)) {
-                                                back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
-                                                                                                                                                                                        (*ri).mMinT,
-                                                                                                                                                                                        (*ri).mRay->mT)
-                                                                                                 );
-                                                front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
-                                                                                                                                                                                         (*ri).mRay->mT,
-                                                                                                                                                                                         (*ri).mMaxT));
-                                        } else {
-                                                back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
-                                                                                                                                                                                        (*ri).mRay->mT,
-                                                                                                                                                                                        (*ri).mMaxT));
-                                                front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
-                                                                                                                                                                                         (*ri).mMinT,
-                                                                                                                                                                                         (*ri).mRay->mT));
-                                        }
-                                } else
-                                        if (side == 1) 
-                                                front->AddRay(*ri);
-                                        else
-                                                back->AddRay(*ri);
+                // first unref ray from the former leaf
+                (*ri).mRay->Unref();
+
+                // Debug << "computed t: " << (*ri).mRay->mT << endl;
+                // determine the side of this ray with respect to the plane
+                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
+                
+                if (side == 0) {
+                  if ((*ri).mRay->HasPosDir(axis)) {
+                        back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
+                                                                                          (*ri).mMinT,
+                                                                                          (*ri).mRay->mT)
+                                                 );
+                        front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
+                                                                                           (*ri).mRay->mT,
+                                                                                           (*ri).mMaxT));
+                  } else {
+                        back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
+                                                                                          (*ri).mRay->mT,
+                                                                                          (*ri).mMaxT));
+                        front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
+                                                                                           (*ri).mMinT,
+                                                                                           (*ri).mRay->mT));
+                  }
+                } else
+                  if (side == 1) 
+                        front->AddRay(*ri);
+                  else
+                        back->AddRay(*ri);
       } else
-                                (*ri).mRay->Unref();
+                (*ri).mRay->Unref();
     }
   } else {
@@ -1020 +1087 @@
     
     for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                                ri != leaf->rays.end();
-                                ri++) {
+                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);
+                // 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();
+                (*ri).mRay->Unref();
     }
   }
-        
+
+  front->SetPvsSize(pvsFront);
+  back->SetPvsSize(pvsBack);
+
   // update stats
   stat.rayRefs -= leaf->rays.size();
@@ -1076 +1146 @@
       //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
       if (in->depth >= minCollapseDepth &&
-          in->lastAccessTime <= maxAccessTime) {
-        released = CollapseSubtree(node, time);
-        break;
+                  in->lastAccessTime <= maxAccessTime) {
+                released = CollapseSubtree(node, time);
+                break;
       }
       
       if (in->back->GetAccessTime() <  
-          in->front->GetAccessTime()) {
-        tstack.push(in->front);
-        tstack.push(in->back);
+                  in->front->GetAccessTime()) {
+                tstack.push(in->front);
+                tstack.push(in->back);
       } else {
-        tstack.push(in->back);
-        tstack.push(in->front);
+                tstack.push(in->back);
+                tstack.push(in->front);
       }
     }
@@ -1106 +1176 @@
 VssTreeNode *
 VssTree::SubdivideLeaf(
-                                                                                         VssTreeLeaf *leaf
-                                                                                         )
+                                           VssTreeLeaf *leaf
+                                           )
 {
   VssTreeNode *node = leaf;
@@ -1113 +1183 @@
   AxisAlignedBox3 leafBBox = GetBBox(leaf);
 
-        static int pass = 0;
-        pass ++;
+  static int pass = 0;
+  pass ++;
         
   // check if we should perform a dynamic subdivision of the leaf
-        if (!TerminationCriteriaSatisfied(leaf)) {
+  if (!TerminationCriteriaSatisfied(leaf)) {
     
-                // memory check and realese...
+        // memory check and realese...
     if (GetMemUsage() > maxTotalMemory) {
       ReleaseMemory( pass );
@@ -1129 +1199 @@
     node =
       SubdivideNode(leaf,
-                                                                                leafBBox,
-                                                                                backBBox,
-                                                                                frontBBox
-                                                                                );
+                                        leafBBox,
+                                        backBBox,
+                                        frontBBox
+                                        );
   }
         
@@ -1142 +1212 @@
 void
 VssTree::UpdateRays(VssRayContainer &remove,
-                                                                                VssRayContainer &add
-                                                                                )
+                                        VssRayContainer &add
+                                        )
 {
   VssTreeLeaf::NewMail();
@@ -1160 +1230 @@
       ri++) {
     if ((*ri)->ScheduledForRemoval())
-//      RemoveRay(*ri, NULL, false);
-// !!! BUG - with true it does not work correctly - aggreated delete
+          //      RemoveRay(*ri, NULL, false);
+          // !!! BUG - with true it does not work correctly - aggreated delete
       RemoveRay(*ri, NULL, true);
     else
@@ -1173 +1243 @@
       ri != add.end();
       ri++) {
-                VssTreeNode::RayInfo info(*ri);
-                if (ClipRay(info, bbox))
-                        AddRay(info);
+        VssTreeNode::RayInfo info(*ri);
+        if (ClipRay(info, bbox))
+          AddRay(info);
   }
 }
@@ -1182 +1252 @@
 void
 VssTree::RemoveRay(VssRay *ray,
-                                                                         vector<VssTreeLeaf *> *affectedLeaves,
-                                                                         const bool removeAllScheduledRays
-                                                                         )
+                                   vector<VssTreeLeaf *> *affectedLeaves,
+                                   const bool removeAllScheduledRays
+                                   )
 {
         
@@ -1211 +1281 @@
       
       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;
+                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--;
-            }
-          }
-        }
+                          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;
-          }
-        }
+                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;
+                  }
+                }
       
     }
@@ -1292 +1362 @@
 void
 VssTree::TraverseInternalNode(
-                                                                                                                        RayTraversalData &data,
-                                                                                                                        stack<RayTraversalData> &tstack)
+                                                          RayTraversalData &data,
+                                                          stack<RayTraversalData> &tstack)
 {
   VssTreeInterior *in =  (VssTreeInterior *) data.node;
@@ -1301 +1371 @@
     // determine the side of this ray with respect to the plane
     int side = in->ComputeRayIntersection(data.rayData,
-                                                                                                                                                                        data.rayData.mRay->mT);
+                                                                                  data.rayData.mRay->mT);
     
   
     if (side == 0) {
       if (data.rayData.mRay->HasPosDir(in->axis)) {
-                                tstack.push(RayTraversalData(in->back,
-                                                                                                                                                 VssTreeNode::RayInfo(data.rayData.mRay,
-                                                                                                                                                                                                                                        data.rayData.mMinT,
-                                                                                                                                                                                                                                        data.rayData.mRay->mT))
-                                                                                );
-                                
-                                tstack.push(RayTraversalData(in->front,
-                                                                                                                                                 VssTreeNode::RayInfo(data.rayData.mRay,
-                                                                                                                                                                                                                                        data.rayData.mRay->mT,
-                                                                                                                                                                                                                                        data.rayData.mMaxT
-                                                                                                                                                                                                                                        ))
-                                                                                );
+                tstack.push(RayTraversalData(in->back,
+                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                  data.rayData.mMinT,
+                                                                                                                  data.rayData.mRay->mT))
+                                        );
+                                
+                tstack.push(RayTraversalData(in->front,
+                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                  data.rayData.mRay->mT,
+                                                                                                                  data.rayData.mMaxT
+                                                                                                                  ))
+                                        );
         
       } else {
-                                tstack.push(RayTraversalData(in->back,
-                                                                                                                                                 VssTreeNode::RayInfo(data.rayData.mRay,
-                                                                                                                                                                                                                                        data.rayData.mRay->mT,
-                                                                                                                                                                                                                                        data.rayData.mMaxT
-                                                                                                                                                                                                                                        ))
-                                                                                );
-                                
-                                tstack.push(RayTraversalData(in->front,
-                                                                                                                                                 VssTreeNode::RayInfo(data.rayData.mRay,
-                                                                                                                                                                                                                                        data.rayData.mMinT,
-                                                                                                                                                                                                                                        data.rayData.mRay->mT))
-                                                                                );
+                tstack.push(RayTraversalData(in->back,
+                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                  data.rayData.mRay->mT,
+                                                                                                                  data.rayData.mMaxT
+                                                                                                                  ))
+                                        );
+                                
+                tstack.push(RayTraversalData(in->front,
+                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
+                                                                                                                  data.rayData.mMinT,
+                                                                                                                  data.rayData.mRay->mT))
+                                        );
                                 
                                 
@@ -1337 +1407 @@
     } else
       if (side == 1)
-                                tstack.push(RayTraversalData(in->front, data.rayData));
+                tstack.push(RayTraversalData(in->front, data.rayData));
       else
-                                tstack.push(RayTraversalData(in->back, data.rayData));
+                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));
+        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));
   }
 }
@@ -1367 +1437 @@
 #endif
 
-//    tstat.collapsedSubtrees++;
-//    tstat.collapseDepths += (int)sroot->depth;
-//    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
+  //    tstat.collapsedSubtrees++;
+  //    tstat.collapseDepths += (int)sroot->depth;
+  //    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
   
   tstack.push(sroot);
@@ -1382 +1452 @@
       VssTreeLeaf *leaf = (VssTreeLeaf *) node;
       for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                                        ri != leaf->rays.end();
-                                        ri++) {
-                                
-                                totalRayCount++;
-                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
-                                        (*ri).mRay->Mail();
-                                        rayCount++;
-                                }
+                  ri != leaf->rays.end();
+                  ri++) {
+                                
+                totalRayCount++;
+                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
+                  (*ri).mRay->Mail();
+                  rayCount++;
+                }
       }
     } else {
@@ -1416 +1486 @@
       
       for(VssTreeNode::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();
+                  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();
                                 
       }
@@ -1440 +1510 @@
   delete sroot;
   
-//   for(VssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
-//       ri != newLeaf->rays.end();
-//       ri++)
-//     (*ri).ray->UnMail(2);
+  //   for(VssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
+  //       ri != newLeaf->rays.end();
+  //       ri++)
+  //     (*ri).ray->UnMail(2);
 
 
@@ -1460 +1530 @@
 #endif
 
-        //  tstat.collapsedNodes += collapsedNodes;
-        //  tstat.collapsedRays += totalRayCount - rayCount;
+  //  tstat.collapsedNodes += collapsedNodes;
+  //  tstat.collapsedRays += totalRayCount - rayCount;
     
   return totalRayCount - rayCount;
@@ -1468 +1538 @@
 
 int
-VssTree::GetPvsSize(VssTreeNode *node, const AxisAlignedBox3 &box) const
-{
-        stack<VssTreeNode *> tstack;
+VssTree::GetPvsSize(const AxisAlignedBox3 &box) const
+{
+  stack<VssTreeNode *> tstack;
   tstack.push(root);
 
-        Intersectable::NewMail();
-        int pvsSize = 0;
+  Intersectable::NewMail();
+  int pvsSize = 0;
         
   while (!tstack.empty()) {
@@ -1482 +1552 @@
  
     if (node->IsLeaf()) {
-                        VssTreeLeaf *leaf = (VssTreeLeaf *)node;
-                        for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                                        ri != leaf->rays.end();
-                                        ri++)
-                                if ((*ri).mRay->IsActive()) {
-                                        Intersectable *object;
+          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
+          for(VssTreeNode::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();
-                                        }
+                  object = (*ri).mRay->mOriginObject;
+                  if (object && !object->Mailed()) {
+                        pvsSize++;
+                        object->Mail();
+                  }
 #endif
-                                        object = (*ri).mRay->mTerminationObject;
-                                        if (object && !object->Mailed()) {
-                                                pvsSize++;
-                                                object->Mail();
-                                        }
-                                }
-                } else {
-                        VssTreeInterior *in = (VssTreeInterior *)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);
-                        }
+                  object = (*ri).mRay->mTerminationObject;
+                  if (object && !object->Mailed()) {
+                        pvsSize++;
+                        object->Mail();
+                  }
                 }
-        }
-        return pvsSize;
+        } else {
+          VssTreeInterior *in = (VssTreeInterior *)node;
+          if (in->axis < 3) {
+                if (box.Max(in->axis) >= in->position )
+                  tstack.push(in->front);
+                                
+                if (box.Min(in->axis) <= in->position )
+                  tstack.push(in->back);
+          } else {
+                // both nodes for directional splits
+                tstack.push(in->front);
+                tstack.push(in->back);
+          }
+        }
+  }
+  return pvsSize;
 }
 
 void
 VssTree::GetRayContributionStatistics(
-                                                                                                                                                        float &minRayContribution,
-                                                                                                                                                        float &maxRayContribution,
-                                                                                                                                                        float &avgRayContribution
-                                                                                                                                                        )
-{
-        stack<VssTreeNode *> tstack;
+                                                                          float &minRayContribution,
+                                                                          float &maxRayContribution,
+                                                                          float &avgRayContribution
+                                                                          )
+{
+  stack<VssTreeNode *> tstack;
   tstack.push(root);
         
-        minRayContribution = 1.0f;
-        maxRayContribution = 0.0f;
-        float sumRayContribution = 0.0f;
-        int leaves = 0;
+  minRayContribution = 1.0f;
+  maxRayContribution = 0.0f;
+  float sumRayContribution = 0.0f;
+  int leaves = 0;
 
   while (!tstack.empty()) {
@@ -1539 +1609 @@
                 
     if (node->IsLeaf()) {
-                        leaves++;
-                        VssTreeLeaf *leaf = (VssTreeLeaf *)node;
-                        float c = leaf->GetAvgRayContribution();
-                        if (c > maxRayContribution)
-                                maxRayContribution = c;
-                        if (c < minRayContribution)
-                                minRayContribution = c;
-                        sumRayContribution += c;
+          leaves++;
+          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
+          float c = leaf->GetAvgRayContribution();
+          if (c > maxRayContribution)
+                maxRayContribution = c;
+          if (c < minRayContribution)
+                minRayContribution = c;
+          sumRayContribution += c;
                         
-                } else {
-                        VssTreeInterior *in = (VssTreeInterior *)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;
+        } else {
+          VssTreeInterior *in = (VssTreeInterior *)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;
 }
 
@@ -1564 +1634 @@
 int
 VssTree::GenerateRays(const float ratioPerLeaf,
-                                                                                        SimpleRayContainer &rays)
-{
-        stack<VssTreeNode *> tstack;
+                                          SimpleRayContainer &rays)
+{
+  stack<VssTreeNode *> tstack;
   tstack.push(root);
         
@@ -1574 +1644 @@
                 
     if (node->IsLeaf()) {
-                        VssTreeLeaf *leaf = (VssTreeLeaf *)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));
-                        }
+          VssTreeLeaf *leaf = (VssTreeLeaf *)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 = VssRay::GetDirection(dirVector.x, dirVector.y);
+                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
+                rays.push_back(SimpleRay(origin, direction));
+          }
                         
-                } else {
-                        VssTreeInterior *in = (VssTreeInterior *)node;
-                        // both nodes for directional splits
-                        tstack.push(in->front);
-                        tstack.push(in->back);
-                }
-        }
-
-        return rays.size();
+        } else {
+          VssTreeInterior *in = (VssTreeInterior *)node;
+          // both nodes for directional splits
+          tstack.push(in->front);
+          tstack.push(in->back);
+        }
+  }
+
+  return rays.size();
 }
 
@@ -1601 +1671 @@
 VssTree::CollectLeaves(vector<VssTreeLeaf *> &leaves)
 {
-        stack<VssTreeNode *> tstack;
+  stack<VssTreeNode *> tstack;
   tstack.push(root);
         
@@ -1609 +1679 @@
                 
     if (node->IsLeaf()) {
-                        VssTreeLeaf *leaf = (VssTreeLeaf *)node;
-                        leaves.push_back(leaf);
-                } else {
-                        VssTreeInterior *in = (VssTreeInterior *)node;
-                        // both nodes for directional splits
-                        tstack.push(in->front);
-                        tstack.push(in->back);
-                }
-        }
+          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
+          leaves.push_back(leaf);
+        } else {
+          VssTreeInterior *in = (VssTreeInterior *)node;
+          // both nodes for directional splits
+          tstack.push(in->front);
+          tstack.push(in->back);
+        }
+  }
+}
+
+bool
+VssTree::ValidLeaf(VssTreeLeaf *leaf) const
+{
+  return leaf->rays.size() > termMinRays/4;
+}
+
+void
+GenerateExtendedConvexCombinationWeights(float &w1,
+                                                                                 float &w2,
+                                                                                 float &w3,
+                                                                                 const float overlap
+                                                                                 )
+{
+  w1 = RandomValue(-overlap, 1.0f + overlap);
+  w2 = RandomValue(-overlap, 1.0f + overlap);
+  w3 = RandomValue(-overlap, 1.0f + overlap);
+        
+  float c = 1.0f/(w1 + w2 + w3);
+  w1 *= c;
+  w2 *= c;
+  w3 *= c;
+}
+
+void
+VssTree::GenerateLeafRays(VssTreeLeaf *leaf,
+                                                  const int numberOfRays,
+                                                  SimpleRayContainer &rays)
+{
+  int nrays = leaf->rays.size();
+  for (int i=0; i < numberOfRays; i++) {
+        // pickup 3 random rays
+        int r1 = Random(nrays-1);
+        int r2 = Random(nrays-1);
+        int r3 = Random(nrays-1);
+                
+        Vector3 o1 = leaf->rays[r1].Extrap(RandomValue(leaf->rays[r1].GetMinT(), 
+                                                                                                   leaf->rays[r1].GetMaxT()));
+
+        Vector3 o2 = leaf->rays[r2].Extrap(RandomValue(leaf->rays[r2].GetMinT(), 
+                                                                                                   leaf->rays[r2].GetMaxT()));
+                
+        Vector3 o3 = leaf->rays[r3].Extrap(RandomValue(leaf->rays[r3].GetMinT(), 
+                                                                                                   leaf->rays[r3].GetMaxT()));
+
+        const float overlap = 0.1;
+                
+        Vector3 origin, direction;
+        bool useExtendedConvexCombination = true;
+        if (useExtendedConvexCombination) {
+          float w1, w2, w3;
+          GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
+          origin = w1*o1 + w2*o2 + w3*o3;
+          direction =
+                w1*leaf->rays[r1].mRay->GetDir() +
+                w2*leaf->rays[r2].mRay->GetDir() +
+                w3*leaf->rays[r3].mRay->GetDir();
+        } else {
+          origin = GetBBox(leaf).GetRandomPoint();
+          Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
+          direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
+        }
+        //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
+        rays.push_back(SimpleRay(origin, direction));
+  }
 }
 
 int
 VssTree::GenerateRays(const int numberOfRays,
-                                                                                        const int numberOfLeaves,
-                                                                                        SimpleRayContainer &rays)
-{
-
-        vector<VssTreeLeaf *> leaves;
-        
-        CollectLeaves(leaves);
-        
-        sort(leaves.begin(),
-                         leaves.end(),
-                         GreaterContribution);
+                                          const int numberOfLeaves,
+                                          SimpleRayContainer &rays)
+{
+
+  vector<VssTreeLeaf *> leaves;
+        
+  CollectLeaves(leaves);
+        
+  sort(leaves.begin(),
+           leaves.end(),
+           GreaterContribution);
                          
 
-        float sumContrib = 0.0;
-        int i;
-        for (i=0; i < numberOfLeaves; i++) {
-                float c = leaves[i]->GetAvgRayContribution();
-                sumContrib += c;
-                //              cout<<"ray contrib "<<i<<" : "<<c<<endl;
-        }
-
-        float avgContrib = sumContrib/numberOfLeaves;
-        float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
-        
-        for (i=0; i < numberOfLeaves; i++) {
-                VssTreeLeaf *leaf = leaves[i];
-                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));
-                }
-        }
-        return rays.size();
+  float sumContrib = 0.0;
+  int i;
+  int k = 0;
+  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
+        if (ValidLeaf(leaves[i])) {
+          float c = leaves[i]->GetAvgRayContribution();
+          sumContrib += c;
+          //            cout<<"ray contrib "<<i<<" : "<<c<<endl;
+          k++;
+        }
+                                
+  float avgContrib = sumContrib/numberOfLeaves;
+  float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
+  k = 0;
+  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
+        if (ValidLeaf(leaves[i])) {
+          k++;
+          VssTreeLeaf *leaf = leaves[i];
+          float c = leaf->GetAvgRayContribution();
+          int num = (c*ratioPerLeaf + 0.5);
+          GenerateLeafRays(leaf, num, rays);
+        }
+
+  return rays.size();
 }
 
@@ -1667 +1801 @@
 VssTree::GetAvgPvsSize()
 {
-        stack<VssTreeNode *> tstack;
+  stack<VssTreeNode *> tstack;
   tstack.push(root);
 
-        int sumPvs = 0;
-        int leaves = 0;
+  int sumPvs = 0;
+  int leaves = 0;
   while (!tstack.empty()) {
     VssTreeNode *node = tstack.top();
@@ -1677 +1811 @@
                 
     if (node->IsLeaf()) {
-                        VssTreeLeaf *leaf = (VssTreeLeaf *)node;
-                        // update pvs size
-                        leaf->UpdatePvsSize();
-                        sumPvs += leaf->GetPvsSize();
-                        leaves++;
-                } else {
-                        VssTreeInterior *in = (VssTreeInterior *)node;
-                        // both nodes for directional splits
-                        tstack.push(in->front);
-                        tstack.push(in->back);
-                }
-        }
-
-
-        return sumPvs/(float)leaves;
-}
-
-        
+          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
+          // update pvs size
+          leaf->UpdatePvsSize();
+          sumPvs += leaf->GetPvsSize();
+          leaves++;
+        } else {
+          VssTreeInterior *in = (VssTreeInterior *)node;
+          // both nodes for directional splits
+          tstack.push(in->front);
+          tstack.push(in->back);
+        }
+  }
+
+
+  return sumPvs/(float)leaves;
+}
+
+        

trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.h

@@ -226 +226 @@
                 return mRay->GetOrigin(axis) + GetMaxT()*mRay->GetDir(axis);
         }
+
+        Vector3 Extrap(const float t) const {
+                return mRay->Extrap(t);
+        }
         
 #if USE_FIXEDPOINT_T
@@ -411 +415 @@
   
   RayInfoContainer rays;
+        int mPassingRays;
+        
         bool mValidPvs;
+        
+        
         
   VssTreeLeaf(VssTreeInterior *p,
                                                         const int nRays
-                                                        ):VssTreeNode(p), rays(), mPvsSize(0), mValidPvs(false) {
+                                                        ):VssTreeNode(p), rays(), mPvsSize(0), mPassingRays(0), mValidPvs(false) {
     rays.reserve(nRays);
   }
@@ -431 +439 @@
     rays.push_back(data);
     data.mRay->Ref();
+                if (data.GetRayClass() == RayInfo::PASSING_RAY)
+                        mPassingRays++;
   }
         
@@ -438 +448 @@
         void SetPvsSize(const int s) {
                 mPvsSize = s;
+                mValidPvs = true;
         }
 
@@ -537 +548 @@
                                 leafb->rays.size()*b.bbox.GetVolume();
 #endif
-#if 1
+#if 0
                         return
                                 leafa->GetPvsSize()*a.bbox.GetVolume()
@@ -555 +566 @@
                                 leafb->GetPvsSize()/(leafb->rays.size()+1);
 #endif
-#if 0
+#if 1
                         return
                                 leafa->GetPvsSize()*leafa->rays.size()
@@ -623 +634 @@
 
   // type of the splitting to use fo rthe tree construction
-  enum {ESplitRegular, ESplitHeuristic };
+  enum {ESplitRegular, ESplitHeuristic, ESplitHybrid };
   int splitType;
+
+        bool mSplitUseOnlyDrivingAxis;
+
+        // use ray space subdivision instead of view space subdivision
+        bool mUseRss;
         
   // maximal size of the box on which the refdir splitting can be performed
@@ -723 +739 @@
   }
   
+
         float
-        BestCostRatioHeuristic(
-                                                                                                 VssTreeLeaf *node,
-                                                                                                 int &axis,
-                                                                                                 float &position,
-                                                                                                 int &raysBack,
-                                                                                                 int &raysFront,
-                                                                                                 int &pvsBack,
-                                                                                                 int &pvsFront
-                                                                                                 );
-
-        float
-        BestCostRatioRegular(
-                                                                                         VssTreeLeaf *node,
-                                                                                         int &axis,
-                                                                                         float &position,
-                                                                                         int &raysBack,
-                                                                                         int &raysFront,
-                                                                                         int &pvsBack,
-                                                                                         int &pvsFront
-
-                                                                                         );
+        BestCostRatio(
+                                                                VssTreeLeaf *node,
+                                                                int &axis,
+                                                                float &position,
+                                                                int &raysBack,
+                                                                int &raysFront,
+                                                                int &pvsBack,
+                                                                int &pvsFront
+                                                                );
         
         float
@@ -757 +762 @@
                                                                 );
 
-  AxisAlignedBox3 GetBBox(const VssTreeNode *node) {
+        float
+        EvalCostRatioHeuristic(
+                                                                                                 VssTreeLeaf *node,
+                                                                                                 const int axis,
+                                                                                                 float &position,
+                                                                                                 int &raysBack,
+                                                                                                 int &raysFront,
+                                                                                                 int &pvsBack,
+                                                                                                 int &pvsFront
+                                                                                                 );
+
+        float
+        GetCostRatio(
+                                                         VssTreeLeaf *leaf,
+                                                         const int axis,
+                                                         const float position,
+                                                         const int raysBack,
+                                                         const int raysFront,
+                                                         const int pvsBack,
+                                                         const int pvsFront
+                                                         );
+
+  AxisAlignedBox3 GetBBox(const VssTreeNode *node) const {
     if (node->parent == NULL)
       return bbox;
@@ -775 +802 @@
   }
 
-  AxisAlignedBox3 GetDirBBox(const VssTreeNode *node) {
+  AxisAlignedBox3 GetDirBBox(const VssTreeNode *node) const {
 
     if (node->parent == NULL)
@@ -834 +861 @@
         int
         GetRootPvsSize() const {
-                return GetPvsSize(root, bbox);
+                return GetPvsSize(bbox);
         }
         
         int
-        GetPvsSize(VssTreeNode *node, const AxisAlignedBox3 &box) const;
+        GetPvsSize(const AxisAlignedBox3 &box) const;
 
         void
@@ -874 +901 @@
                                         );
 
-        
+        VssTreeNode *GetRoot() const { return root; }
+
+        bool
+        ValidLeaf(VssTreeLeaf *leaf) const;
+
+        void
+        GenerateLeafRays(VssTreeLeaf *leaf,
+                                                                         const int numberOfRays,
+                                                                         SimpleRayContainer &rays);
+
+
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp

@@ -10 +10 @@
 #include "VssRay.h"
 #include "VspKdTree.h"
+#include "VssTree.h"
 
 X3dExporter::X3dExporter(const string filename):Exporter(filename)
@@ -154 +155 @@
   for (; ri != rays.end(); ri++) {
     Vector3 a = (*ri)->GetOrigin();
-        Vector3 b = (*ri)->mTerminationObject ? (*ri)->GetTermination() : a + 1000 * Normalize(b - a);
+                Vector3 b = (*ri)->mTerminationObject ? (*ri)->GetTermination() : a + 1000 * Normalize(b - a);
         
     stream<<a.x<<" "<<a.y<<" "<<a.z<<" ,";
-        stream<<b.x<<" "<<b.y<<" "<<b.z<<" ,\n";
+                stream<<b.x<<" "<<b.y<<" "<<b.z<<" ,\n";
   }
   
@@ -699 +700 @@
 }
 
+
+void
+X3dExporter::AddBoxToMesh(const AxisAlignedBox3 &box,
+                                                                                                        Mesh *mesh)
+{
+        // add 6 vertices of the box
+        int index = (int)mesh->mVertices.size();
+        
+        for (int i=0; i < 8; i++) {
+                Vector3 v;
+                box.GetVertex(i, v);
+                mesh->mVertices.push_back(v);
+        }
+        
+        mesh->AddFace(new Face(index + 0, index + 1, index + 3, index + 2) );
+        mesh->AddFace(new Face(index + 0, index + 2, index + 6, index + 4) );
+        mesh->AddFace(new Face(index + 4, index + 6, index + 7, index + 5) );
+        
+        mesh->AddFace(new Face(index + 3, index + 1, index + 5, index + 7) );
+        mesh->AddFace(new Face(index + 0, index + 4, index + 5, index + 1) );
+        mesh->AddFace(new Face(index + 2, index + 3, index + 7, index + 6) );
+
+}
+
+bool
+X3dExporter::ExportVssTree(const VssTree &tree)
+{
+  stack<VssTreeNode *> tStack;
+        
+  tStack.push(tree.GetRoot());
+        
+  Mesh *mesh = new Mesh;
+  VssRayContainer rays;
+        
+  while (!tStack.empty()) {
+
+                VssTreeNode *node = tStack.top();
+    tStack.pop();
+
+                        
+    if (!node->IsLeaf()) {
+      VssTreeInterior *interior = (VssTreeInterior *)node;
+      tStack.push(interior->front);
+      tStack.push(interior->back);
+    } else {
+                        VssTreeLeaf *leaf = (VssTreeLeaf *)node;
+                        AxisAlignedBox3 box;
+                        box = tree.GetBBox(leaf);
+                        AddBoxToMesh(box, mesh);
+
+                        if (tree.ValidLeaf(leaf)) {
+                                
+                                Vector3 origin = box.Center();
+                                box = tree.GetDirBBox(leaf);
+                                VssRay *ray;
+                                
+                                const indices[][2] = {{0,0}, {0,1}, {1,1}, {1,0}}; 
+                                MeshInstance dummy(mesh);
+                                for (int i=0; i < 4; i++) {
+                                        //                              Vector3 v = box.GetVertex(indices[i][0], indices[i][1], 0);
+                                        Vector3 v = box.Center();
+                                        
+                                        Vector3 direction = VssRay::GetDirection(v.x, v.y);
+                                        direction.Normalize();
+                                        float k = 100.0f*leaf->GetAvgRayContribution();
+                                        // get 4 corners of the ray directions
+                                        
+                                        ray = new VssRay(origin, origin + (direction*k), NULL, &dummy);
+                                        rays.push_back(ray);
+                                }
+                        }
+                }
+  }
+
+  ExportMesh(mesh);
+        ExportRays(rays);
+        CLEAR_CONTAINER(rays);
+  delete mesh;
+  return true;
+}
 
 bool

trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h

@@ -38 +38 @@
   bool
   ExportKdTree(const KdTree &tree);
-  
+
+        bool
+  ExportVssTree(const VssTree &tree);
+
   bool 
   ExportVspKdTree(const VspKdTree &tree, const int maxPvs);
@@ -108 +111 @@
 
   bool
-  ExportBspTreeRayDensity(const BspTree &tree);  
+  ExportBspTreeRayDensity(const BspTree &tree);
+
+        void
+        AddBoxToMesh(const AxisAlignedBox3 &box,
+                                                         Mesh *mesh);
+                
 };
  

trunk/VUT/GtpVisibilityPreprocessor/src/common.h

@@ -238 +238 @@
 
 
+
 inline Real sqr(Real a)
 {

trunk/VUT/GtpVisibilityPreprocessor/src/default.env

@@ -13 +13 @@
 #;../data/vienna/vienna-plane.x3d
 # filename ../data/vienna/viewcells-25-sel.x3d
-filename ../data/atlanta/atlanta2.x3d
+# filename ../data/atlanta/atlanta2.x3d
 # filename ../data/soda/soda.dat
-# filename ../data/soda/soda5.dat
+ filename ../data/soda/soda5.dat
 }
 
@@ -25 +25 @@
 VssPreprocessor {
         samplesPerPass  100000
-        initialSamples 2000000
+        initialSamples 200000
         vssSamples 1000000
-        vssSamplesPerPass 100000
+        vssSamplesPerPass 50000
         useImportanceSampling true
 }
 
 VssTree {
+        useRss    true
         epsilon         1e-6
 
@@ -38 +39 @@
         minRays         100
         minSize         0.001
-  maxCostRatio  0.98
+  maxCostRatio  1.1
         maxRayContribution 0.05
         
@@ -44 +45 @@
         maxStaticMemory 50
 
-        splitType regular
-#       splitType heuristics
+#       splitType regular
+#       splitType heuristic
+        splitType hybrid
+        splitUseOnlyDrivingAxis true
 
         numberOfEndPointDomains 10000