Changeset 492

Timestamp:

01/03/06 23:33:45 (18 years ago)

Author:

bittner

Message:

Large merge - viewcells seem not functional now

Location:

trunk/VUT/GtpVisibilityPreprocessor/src

Files:

• AxisAlignedBox3.cpp (modified) (4 diffs)
• AxisAlignedBox3.h (modified) (4 diffs)
• Camera.cpp (modified) (4 diffs)
• Camera.h (modified) (3 diffs)
• Environment.cpp (modified) (19 diffs)
• Environment.h (modified) (4 diffs)
• GlRenderer.cpp (added)
• GlRenderer.h (added)
• Halton.cpp (modified) (1 diff)
• Halton.h (modified) (1 diff)
• Intersectable.h (modified) (2 diffs)
• KdTree.cpp (modified) (1 diff)
• Makefile (modified) (5 diffs)
• Mesh.cpp (modified) (7 diffs)
• Mesh.h (modified) (2 diffs)
• Preprocessor.cpp (modified) (5 diffs)
• Preprocessor.h (modified) (4 diffs)
• PreprocessorThread.cpp (added)
• PreprocessorThread.h (added)
• Pvs.h (modified) (4 diffs)
• RssPreprocessor.cpp (modified) (22 diffs)
• RssPreprocessor.h (modified) (6 diffs)
• RssTree.cpp (modified) (54 diffs)
• RssTree.h (modified) (15 diffs)
• SamplingPreprocessor.cpp (modified) (6 diffs)
• SceneGraph.cpp (modified) (2 diffs)
• SceneGraph.h (modified) (3 diffs)
• Vector3.cpp (modified) (1 diff)
• Vector3.h (modified) (2 diffs)
• ViewCellBsp.cpp (modified) (1 diff)
• ViewCellBsp.h (modified) (1 diff)
• ViewCellsManager.cpp (modified) (4 diffs)
• ViewCellsManager.h (modified) (6 diffs)
• VspBspTree.cpp (modified) (2 diffs)
• VspBspTree.h (modified) (1 diff)
• VssPreprocessor.cpp (modified) (1 diff)
• VssRay.cpp (modified) (3 diffs)
• VssRay.h (modified) (8 diffs)
• X3dExporter.cpp (modified) (1 diff)
• default.env (modified) (9 diffs)
• main.cpp (modified) (4 diffs)
• preprocessor.pro (modified) (3 diffs)

trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp

@@ -301 +301 @@
       maxx = t1;
     }
-    if (maxx < 0.0)
-      return 0;
+        //    if (maxx < 0.0)
+        //      return 0;
   }
 
@@ -327 +327 @@
       maxx = t1;
     }
-    if (maxx < 0.0)
-      return 0;
+        //    if (maxx < 0.0)
+        //      return 0;
   }
 
@@ -355 +355 @@
       maxx = t1;
     }
-    if (maxx < 0.0)
-      return 0;
+        //    if (maxx < 0.0)
+        //      return 0;
   }
 
@@ -1774 +1774 @@
 
 
+  // Compute tmin and tmax for a ray, whenever required .. need not pierce box
+int
+AxisAlignedBox3::ComputeMinMaxT(const Vector3 &origin,
+                                                                const Vector3 &direction,
+                                                                float *tmin,
+                                                                float *tmax) const
+{
+
+  register float minx, maxx;
+
+  
+  Vector3 invDirection;
+  const float eps = 1e-6;
+  const float invEps = 1e6;
+  
+  // it does change the ray direction very slightly,
+  // but the size direction vector is not practically changed
+  
+  if (fabs(direction.x) < eps) {
+    if (direction.x < 0.0)
+      invDirection.x = -invEps;
+    else
+      invDirection.x = invEps;
+  }
+  else
+    invDirection.x = 1.0 / direction.x;
+  
+  if (fabs(direction.y) < eps) {
+    if (direction.y < 0.0)
+      invDirection.y = -invEps;
+    else
+      invDirection.y = invEps;
+  }
+  else
+    invDirection.y = 1.0 / direction.y;
+  
+  if (fabs(direction.z) < eps) {
+    if (direction.z < 0.0)
+      invDirection.z = -invEps;
+    else
+      invDirection.z = invEps;
+  }
+  else
+    invDirection.z = 1.0 / direction.z;
+
+
+  
+  if (fabs(direction.x) < 0.001) {
+    if (mMin.x < origin.x && mMax.x > origin.x) {
+      minx = -MAXFLOAT;
+      maxx = MAXFLOAT;
+    }
+    else
+      return 0;
+  }
+  else {
+    float t1 = (mMin.x - origin.x) * invDirection.x;
+    float t2 = (mMax.x - origin.x) * invDirection.x;
+    if (t1 < t2) {
+      minx = t1;
+      maxx = t2;
+    }
+    else {
+      minx = t2;
+      maxx = t1;
+    }
+        //    if (maxx < 0.0)
+        //      return 0;
+  }
+
+  *tmin = minx;
+  *tmax = maxx;
+  
+  if (fabs(direction.y) < 0.001) {
+    if (mMin.y < origin.y && mMax.y > origin.y) {
+      minx = -MAXFLOAT;
+      maxx = MAXFLOAT;
+    }
+    else
+      return 0;
+  }
+  else {
+    float t1 = (mMin.y - origin.y) * invDirection.y;
+    float t2 = (mMax.y - origin.y) * invDirection.y;
+    if (t1 < t2) {
+      minx = t1;
+      maxx = t2;
+    }
+    else {
+      minx = t2;
+      maxx = t1;
+    }
+        //    if (maxx < 0.0)
+        //      return 0;
+  }
+
+  if (minx > *tmin)
+    *tmin = minx;
+  if (maxx < *tmax)
+    *tmax = maxx;
+  
+  if (fabs(direction.z) < 0.001) {
+    if (mMin.z < origin.z && mMax.z > origin.z) {
+      minx = -MAXFLOAT;
+      maxx = MAXFLOAT;
+    }
+    else
+      return 0;
+  }
+  else {
+    float t1 = (mMin.z - origin.z) * invDirection.z;
+    float t2 = (mMax.z - origin.z) * invDirection.z;
+    if (t1 < t2) {
+      minx = t1;
+      maxx = t2;
+    }
+    else {
+      minx = t2;
+      maxx = t1;
+    }
+        //    if (maxx < 0.0)
+        //      return 0;
+  }
+
+  if (minx > *tmin)
+    *tmin = minx;
+  if (maxx < *tmax)
+    *tmax = maxx;
+
+  return 1; // yes, intersection was found
+}
+
+
 bool AxisAlignedBox3::GetIntersectionFace(Rectangle3 &face, 
                                                                                   const AxisAlignedBox3 &neighbour) const

trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h

@@ -68 +68 @@
   }
 
+
   void SetMin(const Vector3 &v) {
     mMin = v;
@@ -142 +143 @@
   bool IsPoint() const;
 
+  void Scale(const float scale) {
+        Vector3 newSize = Size()*(scale*0.5f);
+        Vector3 center = Center();
+        mMin = center - newSize;
+        mMax = center + newSize;
+  }
+  
   void
   GetSqrDistances(const Vector3 &point,
@@ -181 +189 @@
                                                                                                         RandomValue(0, size.z)); 
   }
-  
+
+
+  Vector3 GetPoint(const Vector3 &p) const {
+    return mMin + p*Size(); 
+  }
+
   // Returns the smallest axis-aligned box that includes all points
   // inside the two given boxes.
@@ -208 +221 @@
                 ID_Right = 4, ID_Top = 5};
   
+  int
+  ComputeMinMaxT(const Vector3 &origin,
+                                 const Vector3 &direction,
+                                 float *tmin,
+                                 float *tmax) const;
+        
   // Compute tmin and tmax for a ray, whenever required .. need not pierce box
   int ComputeMinMaxT(const Ray &ray, float *tmin, float *tmax) const;
 
   // Compute tmin and tmax for a ray, whenever required .. need not pierce box
-  int ComputeMinMaxT(const Ray &ray, float *tmin, float *tmax,
-                     EFaces &entryFace, EFaces &exitFace) const;
+  int ComputeMinMaxT(const Ray &ray,
+                                         float *tmin,
+                                         float *tmax,
+                                         EFaces &entryFace,
+                                         EFaces &exitFace) const;
   
   // If a ray pierces the box .. returns 1, otherwise 0.

trunk/VUT/GtpVisibilityPreprocessor/src/Camera.cpp

@@ -11 +11 @@
 #include "Mesh.h"
 #include "Exporter.h"
+#include "SceneGraph.h"
 
 
@@ -33 +34 @@
 bool
 Camera::SnapImage(string filename,
-                  KdTree *tree
-                  )
+                                  KdTree *tree,
+                                  SceneGraph *sceneGraph
+                                  )
 {
   int x;
@@ -59 +61 @@
       SetupRay(ray, x, mHeight - (y + 1));
       if (tree->CastRay(ray))
-        if (ray.intersections.size()) {
-          sort(ray.intersections.begin(), ray.intersections.end());
-          MeshInstance *mesh = (MeshInstance*)ray.intersections[0].mObject;
-          RgbColor color(1,1,1);
-          if (mesh->GetMesh()->mMaterial)
-            color = mesh->GetMesh()->mMaterial->mDiffuseColor;
-
-          pbuffer[0] = color.r;
-          pbuffer[1] = color.g;
-          pbuffer[2] = color.b;
-          pbuffer[3] = 1.0f;
-
-        }
+                if (ray.intersections.size()) {
+                  sort(ray.intersections.begin(), ray.intersections.end());
+                  MeshInstance *mesh = (MeshInstance*)ray.intersections[0].mObject;
+                  RgbColor color(1,1,1);
+                  if (mesh->GetMesh()->mMaterial)
+                        color = mesh->GetMesh()->mMaterial->mDiffuseColor;
+                  
+                  pbuffer[0] = color.r;
+                  pbuffer[1] = color.g;
+                  pbuffer[2] = color.b;
+                  pbuffer[3] = 1.0f;
+                  
+                }
       pbuffer+=components;
       
-      if (exportRays && (x==222) && (y==97)) {
-                                Ray *nray = new Ray(ray);
-                                rays.push_back(nray);
-                        }
+          //      if (exportRays && ( 1||(x==222) && (y==97))) {
+          //      if (exportRays && ((x%4==0) && (y%4==0))) {
+          if (exportRays && ray.intersections.empty()) {
+                Ray *nray = new Ray(ray);
+                rays.push_back(nray);
+          }
     }
     //    pbuffer-=2*components*mWidth;
@@ -85 +89 @@
   cout<<"#RAY_CAST_TIME\n";
   cout<<TimeDiff(t1, t2)<<"\n";
-
+  
   
   cout<<"Saving image"<<endl;
-
+  
   ilRegisterType(IL_FLOAT);
   ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
   ilSaveImage((char *const)filename.c_str());
   delete buffer;
+
+  cout<<"done."<<endl<<flush;
   
   Exporter *exporter = NULL;
-  if (exportRays) {
+  if (1) {
     exporter = Exporter::GetExporter(filename + "-rays" + ".x3d");
-    exporter->SetWireframe();
-    exporter->ExportKdTree(*tree);
-                //exporter->ExportBspTree(*bsptree);
-    exporter->ExportRays(rays, 10000);
+    exporter->SetFilled();
+        if (sceneGraph)
+          exporter->ExportScene(sceneGraph->mRoot);
+        
+    //exporter->ExportKdTree(*tree);
+        //exporter->ExportBspTree(*bsptree);
+    exporter->ExportRays(rays, 2000);
     exporter->SetFilled();
     int k =0;
     for (int j=0; j < rays.size(); j++)
       if (rays[j]->kdLeaves.size()) {
-                                Ray *ray = rays[j];
-                                int i;
-                                if (1)
-                                        for (i= 0; i < ray->kdLeaves.size(); i++)
-                                                exporter->ExportBox(tree->GetBox(ray->kdLeaves[i]));
-                                if (0)
-                                        for (i= 0; i < ray->testedObjects.size(); i++)
-                                                exporter->ExportIntersectable(ray->testedObjects[i]);
-                                
+                Ray *ray = rays[j];
+                int i;
+                if (0)
+                  for (i= 0; i < ray->kdLeaves.size(); i++)
+                        exporter->ExportBox(tree->GetBox(ray->kdLeaves[i]));
+                if (0)
+                  for (i= 0; i < ray->testedObjects.size(); i++)
+                        exporter->ExportIntersectable(ray->testedObjects[i]);
+                
       }
-                
+        
     delete exporter;
   }
-  
   
   return true;

trunk/VUT/GtpVisibilityPreprocessor/src/Camera.h

@@ -6 +6 @@
 
 class KdTree;
+class SceneGraph;
 
 class Camera
@@ -22 +23 @@
   
   Camera() {
-    mWidth = 640;
-    mHeight = 480;
+    mWidth = 1400;
+    mHeight = 1000;
     mFovy = 60.0f*(float)M_PI/180.0f;
   }
@@ -63 +64 @@
   
   bool
-  SnapImage(string filename, KdTree *tree);
+  SnapImage(string filename,
+                        KdTree *tree,
+                        SceneGraph *sceneGraph
+                        );
 
   void SetupRay(Ray &ray, const int x, const int y);

trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp

@@ -49 +49 @@
 bool
 Environment::CheckForSwitch(const int argc,
-                             const char *argv[],
-                             const char swtch) const
+                                                        char **argv,
+                                                        const char swtch) const
 {
   for (int i = 1; i < argc; i++)
@@ -130 +130 @@
 void
 Environment::ReadCmdlineParams(const int argc,
-                                const char *argv[],
-                                const char *optParams)
+                                                           char **argv,
+                                                           const char *optParams)
 {
   int i;
@@ -619 +619 @@
 void
 Environment::ParseCmdline(const int argc,
-                           const char *argv[],
-                           const int index)
+                                                  char **argv,
+                                                  const int index)
 {
   int curIndex = -1;
@@ -1064 +1064 @@
 
   RegisterOption("Scene.filename",
-                 optString,
-                 "-scene_filename=",
-                 "atlanta2.x3d");
+                                 optString,
+                                 "scene_filename=",
+                                 "atlanta2.x3d");
 
   RegisterOption("Unigraphics.meshGrouping",
-                 optInt,
-                 "-unigraphics_mesh_grouping=",
-                 "0");
+                                 optInt,
+                                 "unigraphics_mesh_grouping=",
+                                 "0");
   
 
   RegisterOption("KdTree.Termination.minCost",
-                 optInt,
-                 "-kd_term_min_cost=",
-                 "10");
+                                 optInt,
+                                 "kd_term_min_cost=",
+                                 "10");
   
   RegisterOption("KdTree.Termination.maxDepth",
-                 optInt,
-                 "-kd_term_max_depth=",
-                 "20");
+                                 optInt,
+                                 "kd_term_max_depth=",
+                                 "20");
 
   RegisterOption("KdTree.Termination.maxCostRatio",
-                 optFloat,
-                 "-kd_term_max_cost_ratio=",
-                 "1.5");
+                                 optFloat,
+                                 "kd_term_max_cost_ratio=",
+                                 "1.5");
 
   RegisterOption("KdTree.Termination.ct_div_ci",
-                 optFloat,
-                 "-kd_term_ct_div_ci=",
-                 "1.0");
+                                 optFloat,
+                                 "kd_term_ct_div_ci=",
+                                 "1.0");
 
   RegisterOption("KdTree.splitMethod",
-                 optString,
-                 "-kd_split_method=",
-                 "spatialMedian");
+                                 optString,
+                                 "kd_split_method=",
+                                 "spatialMedian");
 
   RegisterOption("KdTree.splitBorder",
                  optFloat,
-                 "-kd_split_border=",
+                 "kd_split_border=",
                  "0.1");
 
   RegisterOption("KdTree.sahUseFaces",
                  optBool,
-                 "-kd_sah_use_faces=",
+                 "kd_sah_use_faces=",
                  "true");
 
@@ -1112 +1112 @@
   RegisterOption("MeshKdTree.Termination.minCost",
                  optInt,
-                 "-kd_term_min_cost=",
+                 "kd_term_min_cost=",
                  "10");
   
   RegisterOption("MeshKdTree.Termination.maxDepth",
                  optInt,
-                 "-kd_term_max_depth=",
+                 "kd_term_max_depth=",
                  "20");
 
   RegisterOption("MeshKdTree.Termination.maxCostRatio",
                  optFloat,
-                 "-kd_term_max_cost_ratio=",
+                 "kd_term_max_cost_ratio=",
                  "1.5");
 
   RegisterOption("MeshKdTree.Termination.ct_div_ci",
                  optFloat,
-                 "-kd_term_ct_div_ci=",
+                 "kd_term_ct_div_ci=",
                  "1.0");
 
   RegisterOption("MeshKdTree.splitMethod",
                  optString,
-                 "-kd_split_method=",
+                 "kd_split_method=",
                  "spatialMedian");
 
   RegisterOption("MeshKdTree.splitBorder",
                  optFloat,
-                 "-kd_split_border=",
+                 "kd_split_border=",
                  "0.1");
 
   RegisterOption("Sampling.totalSamples",
                  optInt,
-                 "-total_samples=",
+                 "total_samples=",
                  "1000000");
 
   RegisterOption("Sampling.samplesPerPass",
                  optInt,
-                 "-samples_per_pass=",
+                 "samples_per_pass=",
                  "10");
 
   RegisterOption("VssPreprocessor.initialSamples",
-          optInt,
-          "-initial_samples=",
-          "100000");
-
+                                 optInt,
+                                 "initial_samples=",
+                                 "100000");
+  
   RegisterOption("VssPreprocessor.vssSamples",
-          optInt,
-          "-vss_samples=",
-          "1000000");
+                                 optInt,
+                                 "vss_samples=",
+                                 "1000000");
         
   RegisterOption("VssPreprocessor.vssSamplesPerPass",
-          optInt,
-          "-vss_samples_per_pass=",
-          "1000");
-
+                                 optInt,
+                                 "vss_samples_per_pass=",
+                                 "1000");
+  
   RegisterOption("VssPreprocessor.samplesPerPass",
-          optInt,
-          "-samples_per_pass=",                                                         
-          "100000");
+                                 optInt,
+                                 "samples_per_pass=",
+                                 "100000");
+
 
 
   RegisterOption("VssPreprocessor.useImportanceSampling",
-          optBool,
-          "-vss_use_importance=",
-          "true");
+                                 optBool,
+                                 "vss_use_importance=",
+                                 "true");
 
    RegisterOption("VssPreprocessor.loadInitialSamples",
@@ -1194 +1195 @@
   RegisterOption("ViewCells.type",
                  optString,
-                 "-view_cells_type",
+                 "view_cells_type",
                  "bspTree");
 
   RegisterOption("ViewCells.PostProcess.samples",
-          optInt,
-          "-view_cells_postprocess_samples=",
-          "200000");
-
+                                 optInt,
+                                 "view_cells_postprocess_samples=",
+                                 "200000");
+  
   RegisterOption("ViewCells.Visualization.samples",
-          optInt,
-          "-view_cells_visualization_samples=",
-          "20000");
-
+                                 optInt,
+                                 "view_cells_visualization_samples=",
+                                 "20000");
+  
    RegisterOption("ViewCells.loadFromFile",
-                 optBool,
-                 "-view_cells_load_from_file",
-                 "false");
-
+                                  optBool,
+                                  "view_cells_load_from_file",
+                                  "false");
+   
   RegisterOption("ViewCells.maxViewCells",
                  optInt,
-                 "-view_cells_max_view_cells",
+                 "view_cells_max_view_cells=",
                  "0");
 
   
   RegisterOption("ViewCells.maxPvs",
-                 optInt,
-                 "-view_cells_max_pvs",
-                 "300");
+                                 optInt,
+                                 "view_cells_max_pvs=",
+                                 "300");
 
   RegisterOption("ViewCells.PostProcess.minPvsDif",
-                 optInt,
-                 "-view_cells_post_process_min_pvs_dif",
-                 "10");
+                                 optInt,
+                                 "view_cells_post_process_min_pvs_dif=",
+                                 "10");
 
   RegisterOption("ViewCells.PostProcess.maxPvs",
-                 optInt,
-                 "-view_cells_post_process_max_pvs",
-                 "300");
+                                 optInt,
+                                 "view_cells_post_process_max_pvs=",
+                                 "300");
 
   RegisterOption("ViewCells.PostProcess.minPvs",
-                 optString,
-                 "-view_cells_post_process_min_pvs",
-                 "10");
-
+                                 optString,
+                                 "view_cells_post_process_min_pvs",
+                                 "10");
+  
   RegisterOption("ViewCells.filename",
                  optString,
-                 "-view_cells_filename=",
+                 "view_cells_filename=",
                  "atlanta_viewcells_large.x3d");
 
    RegisterOption("ViewCells.height",
                  optFloat,
-                 "-view_cells_height=",
+                 "view_cells_height=",
                  "5.0");
 
@@ -1270 +1271 @@
   RegisterOption("Simulation.objRenderCost",
                  optFloat,
-                 "-simulation_obj_render_cost",
+                 "simulation_obj_render_cost",
                  "1.0");
 
   RegisterOption("Simulation.vcOverhead",
                  optFloat,
-                 "-simulation_vc_overhead",
+                 "simulation_vc_overhead",
                  "0.05");
 
   RegisterOption("Simulation.moveSpeed",
                  optFloat,
-                 "-simulation_moveSpeed",
+                 "simulation_moveSpeed",
                  "1.0");
 
@@ -1293 +1294 @@
   RegisterOption("BspTree.Construction.input",
           optString,
-          "-bsp_construction_input=",
+          "bsp_construction_input=",
           "fromViewCells");
   
   RegisterOption("BspTree.Construction.samples",
           optInt,
-          "-bsp_construction_samples=",
+          "bsp_construction_samples=",
           "100000");
 
   RegisterOption("BspTree.Construction.epsilon",
           optFloat,
-          "-bsp_construction_side_tolerance=",
+          "bsp_construction_side_tolerance=",
           "0.002");
 
   RegisterOption("BspTree.Termination.minPolygons",
                  optInt,
-                 "-bsp_term_min_polygons=",
+                 "bsp_term_min_polygons=",
                  "5");
 
   RegisterOption("BspTree.Termination.minPvs",
                  optInt,
-                 "-bsp_term_min_pvs=",
+                 "bsp_term_min_pvs=",
                  "20");
 
   RegisterOption("BspTree.Termination.minArea",
                  optFloat,
-                 "-bsp_term_min_area=",
+                 "bsp_term_min_area=",
                  "0.001");
   
   RegisterOption("BspTree.Termination.maxRayContribution",
                  optFloat,
-                 "-bsp_term_ray_contribution=",
+                 "bsp_term_ray_contribution=",
                  "0.005");
 
   RegisterOption("BspTree.Termination.minAccRayLenght",
                  optFloat,
-                 "-bsp_term_min_acc_ray_length=",
+                 "bsp_term_min_acc_ray_length=",
                  "50");
 
    RegisterOption("BspTree.Termination.minRays",
                  optInt,
-                 "-bsp_term_min_rays=",
+                 "bsp_term_min_rays=",
                  "-1");
 
    RegisterOption("BspTree.Termination.ct_div_ci",
                  optFloat,
-                 "-bsp_term_ct_div_ci=",
+                 "bsp_term_ct_div_ci=",
                  "0.0");
   
   RegisterOption("BspTree.Termination.maxDepth",
                  optInt,
-                 "-bsp_term_max_depth=",
+                 "bsp_term_max_depth=",
                  "100");
 
   RegisterOption("BspTree.Termination.maxCostRatio",
                  optFloat,
-                 "-bsp_term_axis_aligned_max_cost_ratio=",
+                 "bsp_term_axis_aligned_max_cost_ratio=",
                  "1.5");
 
   RegisterOption("BspTree.Termination.AxisAligned.ct_div_ci",
                  optFloat,
-                 "-bsp_term_axis_aligned_ct_div_ci=",
+                 "bsp_term_axis_aligned_ct_div_ci=",
                  "0.5");
 
   RegisterOption("BspTree.AxisAligned.splitBorder",
                  optFloat,
-                 "-bsp__axis_aligned_split_border=",
+                 "bsp__axis_aligned_split_border=",
                  "0.1");
 
   RegisterOption("BspTree.Termination.AxisAligned.minPolys",
                  optInt,
-                 "-bsp_term_axis_aligned_max_polygons=",
+                 "bsp_term_axis_aligned_max_polygons=",
                  "50");
 
    RegisterOption("BspTree.Termination.AxisAligned.minObjects",
                  optInt,
-                 "-bsp_term_min_objects=",
+                 "bsp_term_min_objects=",
                  "3");
 
   RegisterOption("BspTree.Termination.AxisAligned.minRays",
                  optInt,
-                 "-bsp_term_axis_aligned_min_rays=",
+                 "bsp_term_axis_aligned_min_rays=",
                  "-1");
 
   RegisterOption("BspTree.splitPlaneStrategy",
                  optString,
-                 "-bsp_split_method=",
+                 "bsp_split_method=",
                  "leastSplits");
 
   RegisterOption("BspTree.maxPolyCandidates",
           optInt,
-          "-bsp_max_poly_candidates=",
+          "bsp_max_poly_candidates=",
           "20");
 
   RegisterOption("BspTree.maxRayCandidates",
           optInt,
-          "-bsp_max_plane_candidates=",
+          "bsp_max_plane_candidates=",
           "20");
 
   RegisterOption("BspTree.maxTests",
           optInt,
-          "-bsp_max_tests=",
+          "bsp_max_tests=",
           "5000");
 
   RegisterOption("BspTree.Visualization.exportSplits",
           optBool,
-          "-bsp_visualization.export_splits",
+          "bsp_visualization.export_splits",
           "false");
 
-
-   RegisterOption("BspTree.Factor.verticalSplits", optFloat, "-bsp_factor_vertical=", "1.0");
-   RegisterOption("BspTree.Factor.largestPolyArea", optFloat, "-bsp_factor_largest_poly=", "1.0");
-   RegisterOption("BspTree.Factor.blockedRays", optFloat, "-bsp_factor_blocked=", "1.0");
-   RegisterOption("BspTree.Factor.leastSplits", optFloat, "-bsp_factor_least_splits=", "1.0");
-   RegisterOption("BspTree.Factor.balancedPolys", optFloat, "-bsp_factor_balanced_polys=", "1.0");
-   RegisterOption("BspTree.Factor.balancedViewCells", optFloat, "-bsp_factor_balanced_view_cells=", "1.0");
-   RegisterOption("BspTree.Factor.leastRaySplits", optFloat, "-bsp_factor_least_ray_splits=", "1.0");
-   RegisterOption("BspTree.Factor.balancedRays", optFloat, "-bsp_factor_balanced_rays=", "1.0");
-   RegisterOption("BspTree.Factor.pvs", optFloat, "-bsp_factor_pvs=", "1.0");
+  RegisterOption("BspTree.Factor.verticalSplits", optFloat, "bsp_factor_vertical=", "1.0");
+   RegisterOption("BspTree.Factor.largestPolyArea", optFloat, "bsp_factor_largest_poly=", "1.0");
+   RegisterOption("BspTree.Factor.blockedRays", optFloat, "bsp_factor_blocked=", "1.0");
+   RegisterOption("BspTree.Factor.leastSplits", optFloat, "bsp_factor_least_splits=", "1.0");
+   RegisterOption("BspTree.Factor.balancedPolys", optFloat, "bsp_factor_balanced_polys=", "1.0");
+   RegisterOption("BspTree.Factor.balancedViewCells", optFloat, "bsp_factor_balanced_view_cells=", "1.0");
+   RegisterOption("BspTree.Factor.leastRaySplits", optFloat, "bsp_factor_least_ray_splits=", "1.0");
+   RegisterOption("BspTree.Factor.balancedRays", optFloat, "bsp_factor_balanced_rays=", "1.0");
+   RegisterOption("BspTree.Factor.pvs", optFloat, "bsp_factor_pvs=", "1.0");
 
    /************************************************************************************/
@@ -1418 +1418 @@
    RegisterOption("Preprocessor.type",
                                  optString,
-                                 "-preprocessor=",
+                                 "preprocessor=",
                                  "sampling");
 
@@ -1426 +1426 @@
                                   "rays.out");
 
-   /**************************************************************************************/
-   /*                  View space partition KD tree related options                      */
-   /**************************************************************************************/
+  RegisterOption("Preprocessor.useGlRenderer",
+                                 optBool,
+                                 "useGlRenderer",
+                                 "false");
+  
+  /**************************************************************************************/
+  /*                  View space partition KD tree related options                      */
+  /**************************************************************************************/
 
    RegisterOption("VspKdTree.Construction.samples",
-           optInt,
-           "-vsp_kd_construction_samples=",
-           "100000");
+                                  optInt,
+                                  "vsp_kd_construction_samples=",
+                                  "100000");
 
    RegisterOption("VspKdTree.Termination.maxDepth", 
@@ -1525 +1530 @@
         RegisterOption("VspKdTree.PostProcess.maxCostRatio",
                         optFloat,
-                        "-vsp_kd_post_process_max_cost_ratio=",
+                        "vsp_kd_post_process_max_cost_ratio=",
                         "0.9");
 
@@ -1568 +1573 @@
         RegisterOption("VssTree.queryType", optString, "qtype=", "static");
 
+
         RegisterOption("VssTree.queryPosWeight", optFloat, "qposweight=", "0.0");
         RegisterOption("VssTree.useRefDirSplits", optBool, "refdir", "false");
@@ -1581 +1587 @@
         RegisterOption("RssPreprocessor.initialSamples",
                                                                         optInt,
-                                                                        "-initial_samples=",
+                                                                        "initial_samples=",
                                                                         "100000");
 
         RegisterOption("RssPreprocessor.vssSamples",
-                                                                        optInt,
-                                                                        "-vss_samples=",
-                                                                        "1000000");
-
+                                   optInt,
+                                   "rss_vss_samples=",
+                                   "1000000");
+        
         RegisterOption("RssPreprocessor.vssSamplesPerPass",
-                                                                        optInt,
-                                                                        "-vss_samples_per_pass=",
-                                                                        "1000");
-
-                RegisterOption("RssPreprocessor.samplesPerPass",
-                                                                        optInt,
-                                                                        "-samples_per_pass=",
-                                                                        "100000");
-
+                                   optInt,
+                                   "rss_vss_samples_per_pass=",
+                                   "1000");
+        
+        RegisterOption("RssPreprocessor.samplesPerPass",
+                                   optInt,
+                                   "rss_samples_per_pass=",
+                                   "100000");
+        
         RegisterOption("RssPreprocessor.useImportanceSampling",
-                                                                        optBool,
-                                                                        "-vss_use_importance=",
-                                                                        "true");
-
+                                   optBool,
+                                   "rss_use_importance",
+                                   "true");
+        
+        RegisterOption("RssPreprocessor.objectBasedSampling",
+                                   optBool,
+                                   "rss_object_based_sampling",
+                                   "true");
+        
+        RegisterOption("RssPreprocessor.directionalSampling",
+                                   optBool,
+                                   "rss_directional_sampling",
+                                   "false");
+        
         RegisterOption("RssTree.maxDepth", optInt, "kd_depth=", "12");
         RegisterOption("RssTree.minPvs", optInt, "kd_minpvs=", "1");
@@ -1609 +1625 @@
         RegisterOption("RssTree.maxCostRatio", optFloat, "maxcost=", "0.95");
         RegisterOption("RssTree.maxRayContribution", optFloat, "maxraycontrib=", "0.5");
-
+        
         RegisterOption("RssTree.epsilon", optFloat, "kd_eps=", "1e-6");
         RegisterOption("RssTree.ct_div_ci", optFloat, "kd_ctdivci=", "1.0");
@@ -1615 +1631 @@
         RegisterOption("RssTree.splitType", optString, "split=", "queries");
         RegisterOption("RssTree.splitUseOnlyDrivingAxis", optBool, "splitdriving=", "false");
-
+        
         RegisterOption("RssTree.numberOfEndPointDomains", optInt, "endpoints=", "10000");
-
+        
         RegisterOption("RssTree.minSize", optFloat, "minsize=", "0.001");
-
+        
         RegisterOption("RssTree.maxTotalMemory", optFloat, "mem=", "60.0");
         RegisterOption("RssTree.maxStaticMemory", optFloat, "statmem=", "8.0");
-
+        
         RegisterOption("RssTree.queryType", optString, "qtype=", "static");
-
+        
         RegisterOption("RssTree.queryPosWeight", optFloat, "qposweight=", "0.0");
         RegisterOption("RssTree.useRefDirSplits", optBool, "refdir", "false");
@@ -1631 +1647 @@
         RegisterOption("RssTree.accessTimeThreshold", optInt, "accesstime=", "1000");
         RegisterOption("RssTree.minCollapseDepth", optInt, "colldepth=", "4");
-
+        
         RegisterOption("RssTree.interleaveDirSplits", optBool, "interleavedirsplits", "true");
         RegisterOption("RssTree.dirSplitDepth", optInt, "dirsplidepth=", "10");
-
+        RegisterOption("RssTree.importanceBasedCost", optBool, "importance_based_cost", "true");
+        RegisterOption("RssTree.maxRays", optInt, "rss_max_rays=", "2000000");
+        
         RegisterOption("RssPreprocessor.Export.pvs", optBool, "rss_export_pvs", "false");
         RegisterOption("RssPreprocessor.Export.rssTree", optBool, "rss_export_rss_tree", "false");
@@ -1641 +1659 @@
         RegisterOption("RssPreprocessor.useViewcells", optBool, "rss_use_viewcells", "false");
         RegisterOption("RssPreprocessor.updateSubdivision", optBool, "rss_update_subdivision", "false");
+        
+        
+  /************************************************************************************/
+  /*               View space partition BSP tree related options                      */
+  /************************************************************************************/
+
 
         RegisterOption("RssPreprocessor.loadInitialSamples",
-          optBool,
-          "-vss_load_loadInitialSamples=",
-          "false");
+                                   optBool,
+                                   "vss_load_loadInitialSamples=",
+                                   "false");
 
         RegisterOption("RssPreprocessor.storeInitialSamples",
-          optBool,
-          "-vss_store_storeInitialSamples=",
-          "false");
+                                   optBool,
+                                   "vss_store_storeInitialSamples=",
+                                   "false");
 
         /************************************************************************************/
         /*               View space partition BSP tree related options                      */
         /************************************************************************************/
-
         RegisterOption("VspBspTree.Termination.minPolygons",
-                optInt,
-                "-vsp_bsp_term_min_polygons=",
-                "5");
+                                   optInt,
+                                   "vsp_bsp_term_min_polygons=",
+                                   "5");
 
         RegisterOption("VspBspTree.Termination.minPvs",
-                optInt,
-                "-vsp_bsp_term_min_pvs=",
-                "20");
-
+                                   optInt,
+                                   "vsp_bsp_term_min_pvs=",
+                                   "20");
         RegisterOption("VspBspTree.Termination.minArea",
-                optFloat,
-                "-vsp_bsp_term_min_area=",
-                "0.001");
+                                   optFloat,
+                                   "vsp_bsp_term_min_area=",
+                                   "0.001");
 
         RegisterOption("VspBspTree.Termination.maxRayContribution",
-                optFloat,
-                "-vsp_bsp_term_ray_contribution=",
-                "0.005");
+                                   optFloat,
+                                   "vsp_bsp_term_ray_contribution=",
+                                   "0.005");
 
         RegisterOption("VspBspTree.Termination.minAccRayLenght",
-                optFloat,
-                "-vsp_bsp_term_min_acc_ray_length=",
-                "50");
+                                   optFloat,
+                                   "vsp_bsp_term_min_acc_ray_length=",
+                                   "50");
 
         RegisterOption("VspBspTree.Termination.minRays",
-                optInt,
-                "-vsp_bsp_term_min_rays=",
-                "-1");
+                                   optInt,
+                                   "vsp_bsp_term_min_rays=",
+                                   "-1");
 
         RegisterOption("VspBspTree.Termination.ct_div_ci",
-                optFloat,
-                "-vsp_bsp_term_ct_div_ci=",
-                "0.0");
+                                   optFloat,
+                                   "vsp_bsp_term_ct_div_ci=",
+                                   "0.0");
 
         RegisterOption("VspBspTree.Termination.maxDepth",
-                optInt,
-                "-vsp_bsp_term_max_depth=",
-                "100");
+                                   optInt,
+                                   "vsp_bsp_term_max_depth=",
+                                   "100");
 
         RegisterOption("VspBspTree.Termination.AxisAligned.maxCostRatio",
@@ -1712 +1734 @@
 
         RegisterOption("VspBspTree.Termination.missTolerance",
-                 optInt,
-                 "-vsp_bsp_term_miss_tolerance=",
-                 "4");
-
+                                   optInt,
+                                   "vsp_bsp_term_miss_tolerance=",
+                                   "4");
         RegisterOption("VspBspTree.splitPlaneStrategy",
-                        optString,
-                        "-vsp_bsp_split_method=",
-                        "leastSplits");
+                                   optString,
+                                   "vsp_bsp_split_method=",
+                                   "leastSplits");
 
         RegisterOption("VspBspTree.maxPolyCandidates",
-                optInt,
-                "-vsp_bsp_max_poly_candidates=",
-                "20");
-
+                                   optInt,
+                                   "vsp_bsp_max_poly_candidates=",
+                                   "20");
         RegisterOption("VspBspTree.maxRayCandidates",
-                optInt,
-                "-vsp_bsp_max_plane_candidates=",
-                "20");
+                                   optInt,
+                                   "vsp_bsp_max_plane_candidates=",
+                                   "20");
 
         RegisterOption("VspBspTree.maxTests",
-                optInt,
-                "-vsp_bsp_max_tests=",
-                "5000");
+                                   optInt,
+                                   "vsp_bsp_max_tests=",
+                                   "5000");
 
         RegisterOption("VspBspTree.Construction.samples",
-                optInt,
-                "-bsp_construction_samples=",
-                "100000");
+                                   optInt,
+                                   "bsp_construction_samples=",
+                                   "100000");
 
         RegisterOption("VspBspTree.Construction.epsilon",
-                optFloat,
-                "-vsp_bsp_construction_side_tolerance=",
-                "0.002");
+                                   optFloat,
+                                   "vsp_bsp_construction_side_tolerance=",
+                                   "0.002");
 
         RegisterOption("VspBspTree.Visualization.exportSplits",
-                optBool,
-                "-vsp_bsp_visualization.export_splits",
-                "false");
+                                   optBool,
+                                   "vsp_bsp_visualization.export_splits",
+                                   "false");
         
         RegisterOption("VspBspTree.splitUseOnlyDrivingAxis", 
-                optBool, 
-                "vsp_bsp_splitdriving=", 
-                "false");
-
-
+                                   optBool, 
+                                   "vsp_bsp_splitdriving=", 
+                                   "false");
+        
         RegisterOption("VspBspTree.Factor.leastRaySplits", optFloat, "-vsp_bsp_factor_least_ray_splits=", "1.0");
         RegisterOption("VspBspTree.Factor.balancedRays", optFloat, "-vsp_bsp_factor_balanced_rays=", "1.0");
@@ -1801 +1820 @@
 
 void
-Environment::Parse(const int argc, const char *argv[], bool useExePath)
+Environment::Parse(const int argc, char **argv, bool useExePath)
 {
   

trunk/VUT/GtpVisibilityPreprocessor/src/Environment.h

@@ -206 +206 @@
    */
   bool CheckForSwitch(const int argc,
-                      const char *argv[],
-                      const char swtch) const;
+                                          char **argv,
+                                          const char swtch) const;
   /** First pass of parsing command line.
     This function parses command line and gets from there all non-optional
@@ -220 +220 @@
    */
   void ReadCmdlineParams(const int argc,
-                         const char *argv[],
-                         const char *optParams);
+                                                 char **argv,
+                                                 const char *optParams);
   /** Reading of the environment file.
     This function reads the environment file.
@@ -236 +236 @@
                  be read in.
    */
-  void ParseCmdline(const int argc, const char *argv[], const int index);
+  void ParseCmdline(const int argc, char **argv, const int index);
 
   /** Clears the environment data structures.
@@ -335 +335 @@
   //@}
 
-  void Parse(const int argc, const char *argv[], bool useExePath);
+  void Parse(const int argc, char **argv, bool useExePath);
 
   void

trunk/VUT/GtpVisibilityPreprocessor/src/Halton.cpp

@@ -5 +5 @@
 
 Halton2 halton2;
+float Halton2::_invBases[2];

trunk/VUT/GtpVisibilityPreprocessor/src/Halton.h

@@ -2 +2 @@
 #define __HALTON_H
 
+#include <iostream>
+using namespace std;
+
 class Halton2 {
-        
-        float _invBases[2];
-        float _prev[2];
+  static float _invBases[2];
+  float _prev[2];
+  
+  float halton(float baseRec, float prev) const {
+        float r = 1 - prev - 1e-10;
+        if (baseRec < r)
+          return prev + baseRec;
+        float h = baseRec;
+        float hh;
+        do {
+          hh = h;
+          h *= baseRec;
+        } while (h >= r);
+        return prev + hh + h - 1;
+  }
+  
+public:
+  
+  void Reset() {
+        _prev[0] =_prev[1] = 0;
+  }
+  
+  Halton2() {
+        _invBases[0] = 1./2;
+        _invBases[1] = 1./3;
+        Reset();
+  }
+  
+  void
+  GetNext(float &a, float &b) {
+        a = halton(_invBases[0], _prev[0]);
+        b = halton(_invBases[1], _prev[1]);
+        _prev[0] = a;
+        _prev[1] = b;
+  }
+};
 
-        float halton(float baseRec, float prev) const {
-                float r = 1 - prev - 1e-10;
-                if (baseRec < r) return prev + baseRec;
-                float h = baseRec;
-                float hh;
-                do {
-                        hh = h;
-                        h *= baseRec;
-                } while (h >= r);
-                return prev + hh + h - 1;
+
+/**
+ * Assert whether the argument is a prime number.
+ * @param number the number to be checked
+ */
+inline bool IsPrime(const int number) {
+  bool isIt = true;
+  for(int i = 2; i < number; i++) {
+        if(number % i == 0) {
+          isIt = false;
+          break;
         }
-        
+  }
+        if(number == 2) {
+          isIt = false;
+        }
+        return isIt;
+}
+
+/**
+ * Find the nth prime number.
+ * @param index the ordinal position in the sequence
+ */
+inline int FindPrime(const int index) {
+  if(index < 1) {
+        cerr<<"FindPrime: The argument must be non-negative."<<endl;
+        return -1;
+  }
+  int prime = 1;
+  int found = 1;
+  while(found != index) {
+        prime += 2;
+          if(IsPrime(prime) == true) {
+                found++;
+          }
+  }
+  return prime;
+}
+
+struct HaltonSequence {
 public:
+  int index;
 
-        void Reset() {
-                _prev[0] =_prev[1] = 0;
+  HaltonSequence():index(1) {}
+
+  void Reset() {
+        index = 1;
+  }
+
+  void GenerateNext() {
+        index++;
+  }
+  
+  /**
+   * Returns the nth number in the sequence, taken from a specified dimension.
+   * @param index the ordinal position in the sequence
+   * @param dimension the dimension
+   */
+
+  double GetNumber(const int dimension)  {
+        int base = FindPrime(dimension);
+        if(base == 1) {
+          base++;  //The first dimension uses base 2.
+        } 
+        double remainder;
+        double output = 0.0;
+        double fraction = 1.0 / (double)base;
+        int N1 = 0;
+        int copyOfIndex = index;
+        if(base >= 2 & index >= 1) {
+          while(copyOfIndex > 0) {
+                N1 = (copyOfIndex / base);
+                remainder = copyOfIndex % base;
+                output += fraction * remainder;
+                copyOfIndex = (int)(copyOfIndex / base);
+                fraction /= (double)base;
+          } 
+          return output;
         }
-
-        Halton2() {
-                _invBases[0] = 1./2;
-                _invBases[1] = 1./3;
-                Reset();
+        else {
+          cerr<<"Error generating Halton sequence."<<endl;
+          exit(1);
         }
-        
-        void
-        GetNext(float &a, float &b) {
-                a = halton(_invBases[0], _prev[0]);
-                b = halton(_invBases[1], _prev[1]);
-                _prev[0] = a;
-                _prev[1] = b;
-        }
+  }
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h

@@ -5 +5 @@
 #include "Pvs.h"
 
+class VssRayContainer;
 
 class Intersectable {
@@ -54 +55 @@
   virtual float IntersectionComplexity() = 0;
   
-        virtual int NumberOfFaces() const = 0;
+  virtual int NumberOfFaces() const = 0;
 
-        virtual int Type() const = 0;
-
+  virtual int Type() const = 0;
+  
   virtual int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal) = 0;
-
-        virtual int
-        GetRandomVisibleSurfacePoint(Vector3 &point,
-                                                                                                                         Vector3 &normal,
-                                                                                                                         const Vector3 &viewpoint,
+  
+  virtual int
+  GetRandomVisibleSurfacePoint(Vector3 &point,
+                                                           Vector3 &normal,
+                                                           const Vector3 &viewpoint,
                                                                                                                          const int maxTries
-                                                                                                                         ) = 0;
-
-        virtual ostream &Describe(ostream &s) = 0;
-        
+                                                           ) = 0;
+  
+  virtual ostream &Describe(ostream &s) = 0;
+  
+  virtual int GenerateSilhouetteRays(const int nrays,
+                                                                         const AxisAlignedBox3 &originBox,
+                                                                         const AxisAlignedBox3 &directionBox,
+                                                                         VssRayContainer &rays
+                                                                         ) {return 0;}
+  
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.cpp

@@ -511 +511 @@
   float maxt = 1e6;
   float mint = 0;
-
+ 
   Intersectable::NewMail();
 

trunk/VUT/GtpVisibilityPreprocessor/src/Makefile

@@ -1 +1 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (2.00a) (Qt 4.1.0-rc1) on: st 14. XII 10:37:40 2005
+# Generated by qmake (2.00a) (Qt 4.1.0-rc1) on: út 3. I 23:13:25 2006
 # Project:  preprocessor.pro
 # Template: app
@@ -71 +71 @@
                 D:\Qt\4.1.0-rc1\mkspecs\features\debug_and_release.prf \
                 D:\Qt\4.1.0-rc1\mkspecs\features\default_post.prf \
+                D:\Qt\4.1.0-rc1\mkspecs\features\qt.prf \
+                D:\Qt\4.1.0-rc1\mkspecs\features\win32\opengl.prf \
+                D:\Qt\4.1.0-rc1\mkspecs\features\moc.prf \
                 D:\Qt\4.1.0-rc1\mkspecs\features\win32\thread.prf \
                 D:\Qt\4.1.0-rc1\mkspecs\features\warn_off.prf \
@@ -78 +81 @@
                 D:\Qt\4.1.0-rc1\mkspecs\features\win32\stl.prf \
                 D:\Qt\4.1.0-rc1\mkspecs\features\shared.prf \
-                D:\Qt\4.1.0-rc1\mkspecs\features\qt.prf \
-                D:\Qt\4.1.0-rc1\mkspecs\features\moc.prf \
                 D:\Qt\4.1.0-rc1\mkspecs\features\resources.prf \
                 D:\Qt\4.1.0-rc1\mkspecs\features\uic.prf
@@ -91 +92 @@
 D:\Qt\4.1.0-rc1\mkspecs\features\debug_and_release.prf:
 D:\Qt\4.1.0-rc1\mkspecs\features\default_post.prf:
+D:\Qt\4.1.0-rc1\mkspecs\features\qt.prf:
+D:\Qt\4.1.0-rc1\mkspecs\features\win32\opengl.prf:
+D:\Qt\4.1.0-rc1\mkspecs\features\moc.prf:
 D:\Qt\4.1.0-rc1\mkspecs\features\win32\thread.prf:
 D:\Qt\4.1.0-rc1\mkspecs\features\warn_off.prf:
@@ -98 +102 @@
 D:\Qt\4.1.0-rc1\mkspecs\features\win32\stl.prf:
 D:\Qt\4.1.0-rc1\mkspecs\features\shared.prf:
-D:\Qt\4.1.0-rc1\mkspecs\features\qt.prf:
-D:\Qt\4.1.0-rc1\mkspecs\features\moc.prf:
 D:\Qt\4.1.0-rc1\mkspecs\features\resources.prf:
 D:\Qt\4.1.0-rc1\mkspecs\features\uic.prf:

trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp

@@ -103 +103 @@
         
   for ( ;
-                                faceIndex < mFaces.size();
-                                faceIndex++) {
+                faceIndex < mFaces.size();
+                faceIndex++) {
     hits += CastRayToFace(faceIndex, ray, nearestT, nearestFace, instance);
     if (mIsConvex && nearestFace != -1)
@@ -122 +122 @@
 int
 Mesh::CastRayToSelectedFaces(
-                                                                                                                 Ray &ray,
-                                                                                                                 const vector<int> &faces,
-                                                                                                                 Intersectable *instance
-                                                                                                                 )
+                                                         Ray &ray,
+                                                         const vector<int> &faces,
+                                                         Intersectable *instance
+                                                         )
 {
   vector<int>::const_iterator fi;
@@ -160 +160 @@
 inline int
 int_lineseg(float px,
-                                                float py,
-                                                float u1,
-                                                float v1,
-                                                float u2,
-                                                float v2)
-{
-  float t;
+                        float py,
+                        float u1,
+                        float v1,
+                        float u2,
+                        float v2)
+{
   float ydiff;
 
@@ -183 +182 @@
   if (fabs(ydiff) < Limits::Small) {      // denominator near 0
     if (((fabs(v1) > Limits::Small) ||
-                                 (u1 > 0) || (u2 > 0)))
+                 (u1 > 0) || (u2 > 0)))
       return 0;
     return 1;
   }
-
-  t = -v1 / ydiff;                // Compute parameter
-
-  return (u1 + t * (u2 - u1)) > 0;
+  
+  double t = -v1 / ydiff;                 // Compute parameter
+
+  double thresh;
+  if (ydiff < 0.0f)
+        thresh = -1e20;
+  else
+        thresh = 1e20;
+
+  
+  return (u1 + t * (u2 - u1)) > thresh; //-Limits::Small;
 }
 
@@ -198 +204 @@
 int
 Mesh::RayFaceIntersection(const int faceIndex,
-                                                                                                        const Ray &ray,
-                                                                                                        float &t,
-                                                                                                        const float nearestT
-                                                                                                        )
+                                                  const Ray &ray,
+                                                  float &t,
+                                                  const float nearestT
+                                                  )
 {
   Face *face  = mFaces[faceIndex];
@@ -244 +250 @@
       mVertices[face->mVertexIndices[i]].ExtractVerts(&u2, &v2, paxis);
       // line u1, v1, u2, v2
-      if ((v2 - v1)*(u1 - u) > (u2 - u1)*(v1 - v))
-                                return Ray::NO_INTERSECTION;
+      if ((v1 - v2)*(u - u1) + (u2 - u1)*(v - v1) > 0)
+                return Ray::NO_INTERSECTION;
       u1 = u2;
       v1 = v2;
@@ -416 +422 @@
 int
 TransformedMeshInstance::CastRay(
-                                                                                                                                 Ray &ray
-                                                                                                                                 )
+                                                                 Ray &ray
+                                                                 )
 {
   ray.ApplyTransform(Invert(mWorldTransform));

trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h

@@ -145 +145 @@
   int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
 
+  int
+  GetRandomVisibleSurfacePoint(Vector3 &point,
+                                                           Vector3 &normal,
+                                                           const Vector3 &viewpoint,
+                                                           const int maxTries
+                                                           );
+  
+  virtual ostream &Describe(ostream &s) {
+        return s<<"Mesh #vertices="<<(int)mVertices.size()<<" #faces="<<(int)mFaces.size();
+  }
+  
+
+};
+
+
+class MeshInstance : public Intersectable {
+protected:
+  Mesh *mMesh;
+  
+public:
+  MeshInstance(Mesh *mesh):Intersectable(), mMesh(mesh)
+  {
+  }
+
+  int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
+
         int
         GetRandomVisibleSurfacePoint(Vector3 &point,
@@ -151 +177 @@
                                                                                                                          const int maxTries
                                                                                                                          );
-                
-        virtual ostream &Describe(ostream &s) {
-                return s<<"Mesh #vertices="<<(int)mVertices.size()<<" #faces="<<(int)mFaces.size();
-        }
-
-};
-
-class MeshInstance : public Intersectable {
-protected:
-  Mesh *mMesh;
-  
-public:
-  MeshInstance(Mesh *mesh):Intersectable(), mMesh(mesh)
-  {
-  }
-
-  int GetRandomSurfacePoint(Vector3 &point, Vector3 &normal);
-
-        int
-        GetRandomVisibleSurfacePoint(Vector3 &point,
-                                                                                                                         Vector3 &normal,
-                                                                                                                         const Vector3 &viewpoint,
-                                                                                                                         const int maxTries
-                                                                                                                         );
 
 

trunk/VUT/GtpVisibilityPreprocessor/src/Preprocessor.cpp

@@ -12 +12 @@
 #include "RenderSimulator.h"
 
+Preprocessor *preprocessor;
+
 Preprocessor::Preprocessor():
 mKdTree(NULL),
@@ -19 +21 @@
 mViewCellsManager(NULL)
 {
+  environment->GetBoolValue("Preprocessor.useGlRenderer", mUseGlRenderer);
 }
 
@@ -55 +58 @@
 
   mSceneGraph = new SceneGraph;
-
   
   Parser *parser;
@@ -95 +97 @@
 
         if (result) {
-                mSceneGraph->AssignObjectIds();
-                int intersectables, faces;
-                mSceneGraph->GetStatistics(intersectables, faces);
-                cout<<filename<<" parsed successfully."<<endl;
-                cout<<"#NUM_OBJECTS (Total numner of objects)\n"<<intersectables<<endl;
-                cout<<"#NUM_FACES (Total numner of faces)\n"<<faces<<endl;
-        }
-
-        
-  return result;
+          
+          mSceneGraph->AssignObjectIds();
+          int intersectables, faces;
+          mSceneGraph->GetStatistics(intersectables, faces);
+          cout<<filename<<" parsed successfully."<<endl;
+          cout<<"#NUM_OBJECTS (Total numner of objects)\n"<<intersectables<<endl;
+          cout<<"#NUM_FACES (Total numner of faces)\n"<<faces<<endl;
+          mSceneGraph->CollectObjects(&mObjects);
+          mSceneGraph->mRoot->UpdateBox();
+        }
+        
+        
+        return result;
 }
 
@@ -139 +144 @@
 bool
 Preprocessor::Export( const string filename,
-                                                                                        const bool scene,
-                                                                                        const bool kdtree,
-                                                                                        const bool bsptree
-                                                                                        )
+                                          const bool scene,
+                                          const bool kdtree,
+                                          const bool bsptree
+                                          )
 {
   Exporter *exporter = Exporter::GetExporter(filename);

trunk/VUT/GtpVisibilityPreprocessor/src/Preprocessor.h

@@ -7 +7 @@
 #include "Mesh.h"
 #include "KdTree.h"
+
+#include <QObject>
 
 class RenderSimulator;
@@ -30 +32 @@
       viewcell loading/generation and the visibility computation itself.
   */
-class Preprocessor {
+class Preprocessor : public QObject {
+  Q_OBJECT
+
 public:
         /** Default constructor initialising e.g., KD tree and BSP tree.
@@ -94 +98 @@
   /// scene graph loaded from file
   SceneGraph *mSceneGraph;
+
+  /// raw array of objects
+  ObjectContainer mObjects;
   
   /// kD-tree organizing the scene graph (occluders + occludees) + viewcells
@@ -115 +122 @@
   VspKdTree *mVspKdTree;
 
+  bool mUseGlRenderer;
+  
+protected:
+
+  /////////////////////////
+
+  /// samples used for construction of the BSP view cells tree.
+  int mBspConstructionSamples;
+   /// samples used for construction of the VSP KD tree.
+  int mVspKdConstructionSamples;
   /** Simulates rendering of the scene.
   */
   RenderSimulator *mRenderSimulator;
+
+  signals:
+  void EvalPvsStat();
+
 };
 
 
 
-
+extern Preprocessor *preprocessor;
 
 #endif

trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h

@@ -60 +60 @@
   */
   PvsData<T> *Find(T sample);
+
+  bool GetSampleContribution(T sample, float &contribution);
   
   /** Adds sample to PVS. 
@@ -71 +73 @@
   */
   int AddSample(T sample);
-
+  
   /** Adds one pvs to another one.
           @returns new pvs size
@@ -157 +159 @@
 
 template <typename T>
-int Pvs<T>::AddSample(T sample)
-{
-  float dummy;
-  return AddSample(sample, dummy) ? 1 : 0;
-}
-
-template <typename T>
-bool Pvs<T>::AddSample(T sample, float &contribution)
+int
+Pvs<T>::AddSample(T sample)
+{
+  PvsData<T> *data = Find(sample);
+  
+  if (data)  {
+        return ++data->mVisibleSamples;
+  }
+  else {
+        mEntries[sample] = PvsData<T>(1);
+        return 1;
+  }
+}
+
+template <typename T>
+bool
+Pvs<T>::AddSample(T sample, float &contribution)
 {
   PvsData<T> *data = Find(sample);
@@ -181 +192 @@
 
 template <typename T>
+bool
+Pvs<T>::GetSampleContribution(T sample, float &contribution) 
+{
+  PvsData<T> *data = Find(sample);
+  
+  if (data)  {
+        contribution = 1.0f/(data->mVisibleSamples + 1);
+        return false;
+  }
+  else {
+        contribution = 1.0f;
+        return true;
+  }
+}
+
+template <typename T>
 bool Pvs<T>::RemoveSample(T sample, const int visibleSamples)
 {

trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp

@@ -11 +11 @@
 #include "ViewCellsManager.h"
 #include "RenderSimulator.h"
+#include "GlRenderer.h"
 
 static bool useViewSpaceBox = false;
@@ -33 +34 @@
   environment->GetIntValue("RssPreprocessor.Export.numRays", mExportNumRays);
   environment->GetBoolValue("RssPreprocessor.useViewcells", mUseViewcells);
-  environment->GetBoolValue("RssPreprocessor.useViewcells", mUseViewcells);
+  environment->GetBoolValue("RssPreprocessor.objectBasedSampling", mObjectBasedSampling);
+  environment->GetBoolValue("RssPreprocessor.directionalSampling", mDirectionalSampling);
 
   environment->GetBoolValue("RssPreprocessor.loadInitialSamples", mLoadInitialSamples);
@@ -39 +41 @@
   
   mStats.open("stats.log");
+
 }
 
 RssPreprocessor::~RssPreprocessor()
 {
-  CLEAR_CONTAINER(mVssRays);
+  // mVssRays get deleted in the tree
+  //  CLEAR_CONTAINER(mVssRays);
 }
 
@@ -78 +82 @@
   Vector3 pointA, pointB;
   float bsize = Magnitude(box.Size());
+
   
   if (mKdTree->CastRay(ray)) {
@@ -86 +91 @@
         // 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);
+        if (box.ComputeMinMaxT(ray, &tmin, &tmax) && tmin < tmax)
+          pointA = ray.Extrap(tmax);
+        else
+          return 0;
   }
 
@@ -108 +111 @@
         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);
-  }
-
+        if (box.ComputeMinMaxT(ray, &tmin, &tmax) && tmin < tmax)
+          pointB = ray.Extrap(tmax);
+        else
+          return 0;
+  }
+
+  //  if (objectA == NULL && objectB != NULL) {
+  if (1) {
+        // cast again to ensure that there is no objectA
+        SetupRay(ray, pointB, direction);
+        if (mKdTree->CastRay(ray)) {
+          objectA = ray.intersections[0].mObject;
+          pointA = ray.Extrap(ray.intersections[0].mT);
+        }
+  }
+  
+  
   VssRay *vssRay  = NULL;
 
-  bool validSample = true;
-  if (detectEmptyViewSpace) {
+  bool validSample = (objectA != objectB);
+  if (0 && detectEmptyViewSpace) {   // consider all samples valid
         // check if the viewpoint lies on the line segment AB
         if (Distance(pointA, pointB) <
@@ -176 +188 @@
   Vector3 point;
   if (!use2dSampling) {
-        Vector3 normal;
-        int i = (int)RandomValue(0, (Real)((int)mObjects.size()-1));
-        Intersectable *object = mObjects[i];
-        object->GetRandomSurfacePoint(point, normal);
+        if (mObjectBasedSampling) {
+          Vector3 normal;
+          int i = RandomValue(0, mObjects.size()-1);
+          Intersectable *object = mObjects[i];
+          object->GetRandomSurfacePoint(point, normal);
+        } else {
+          // select
+          if (mDirectionalSampling) {
+                point = viewpoint + UniformRandomVector();
+          } else {
+                // select the other point uniformly distruted in the whole viewspace
+                AxisAlignedBox3 box;
+                
+                if (viewSpaceBox)
+                  box =*viewSpaceBox;
+                else 
+                  box = mKdTree->GetBox();
+                
+                point = box.GetRandomPoint();
+          }
+        }
   } else {
         AxisAlignedBox3 box;
@@ -187 +216 @@
         else 
           box = mKdTree->GetBox();
-                
+        
         point = box.GetRandomPoint();
         point.y = viewpoint.y;
   }
         
+  return point - viewpoint;
+}
+
+Vector3
+RssPreprocessor::InitialGetDirection(const Vector3 &viewpoint,
+                                                                         AxisAlignedBox3 *viewSpaceBox
+                                                                         )
+{
+  Vector3 point;
+  if (!use2dSampling) {
+        if (1 || mObjectBasedSampling) {
+          Vector3 normal;
+          int i = RandomValue(0, mObjects.size()-1);
+          Intersectable *object = mObjects[i];
+          object->GetRandomSurfacePoint(point, normal);
+        } else {
+          // select
+          if (1 || mDirectionalSampling) {
+                point = viewpoint + UniformRandomVector();
+          } else {
+                // select the other point uniformly distruted in the whole viewspace
+                AxisAlignedBox3 box;
+                
+                if (viewSpaceBox)
+                  box =*viewSpaceBox;
+                else 
+                  box = mKdTree->GetBox();
+                
+                point = box.GetRandomPoint();
+          }
+        }
+  } else {
+        AxisAlignedBox3 box;
+        
+        if (viewSpaceBox)
+          box =*viewSpaceBox;
+        else 
+          box = mKdTree->GetBox();
+        
+        point = box.GetRandomPoint();
+        point.y = viewpoint.y;
+  }
+  
   return point - viewpoint;
 }
@@ -203 +275 @@
   int num;
 
-
-  
-  float avgPvsSize;
-  float avgRayContribution;
-  float avgPvsEntropy;
-  float avgRayLengthEntropy;
-  float avgImportance;
-  float avgRays;
-  rssTree->GetTreeStatistics(
-                                                         avgPvsSize,
-                                                         avgRays,
-                                                         avgRayContribution,
-                                                         avgPvsEntropy,
-                                                         avgRayLengthEntropy,
-                                                         avgImportance);
-  
+  rssTree->UpdateTreeStatistics();
+
   cout<<
-        "#RSS_AVG_PVS_SIZE\n"<<avgPvsSize<<endl<<
-        "#RSS_AVG_RAYS\n"<<avgRays<<endl<<
-        "#RSS_AVG_RAY_CONTRIB\n"<<avgRayContribution<<endl<<
-        "#RSS_AVG_PVS_ENTROPY\n"<<avgPvsEntropy<<endl<<
-        "#RSS_AVG_RAY_LENGTH_ENTROPY\n"<<avgRayLengthEntropy<<endl<<
-        "#RSS_AVG_IMPORTANCE\n"<<avgImportance<<endl;
+        "#RSS_AVG_PVS_SIZE\n"<<rssTree->stat.avgPvsSize<<endl<<
+        "#RSS_AVG_RAYS\n"<<rssTree->stat.avgRays<<endl<<
+        "#RSS_AVG_RAY_CONTRIB\n"<<rssTree->stat.avgRayContribution<<endl<<
+        "#RSS_AVG_PVS_ENTROPY\n"<<rssTree->stat.avgPvsEntropy<<endl<<
+        "#RSS_AVG_RAY_LENGTH_ENTROPY\n"<<rssTree->stat.avgRayLengthEntropy<<endl<<
+        "#RSS_AVG_IMPORTANCE\n"<<rssTree->stat.avgImportance<<endl;
   
   if (0) {
-        float p = desiredSamples/(float)(avgRayContribution*rssTree->stat.Leaves());
+        float p = desiredSamples/(float)(rssTree->stat.avgRayContribution*rssTree->stat.Leaves());
         num = rssTree->GenerateRays(p, rays);
   } else {
@@ -249 +307 @@
   cout<<"Exporting vss rays..."<<endl<<flush;
         
-  float prob = number/(float)vssRays.size();
-
-
   Exporter *exporter = NULL;
   exporter = Exporter::GetExporter(filename);
@@ -266 +321 @@
   }
         
-  VssRayContainer rays; for (int i=0; i < vssRays.size(); i++)
-        if (RandomValue(0,1) < prob)
-          rays.push_back(vssRays[i]);
+  VssRayContainer rays;
+  
+  vssRays.SelectRays( number, rays);
 
   exporter->ExportRays(rays, RgbColor(1, 0, 0));
@@ -388 +443 @@
 }
 
+
+void
+RssPreprocessor::ComputeRenderError()
+{
+  // compute rendering error
+  if (renderer) {
+        emit EvalPvsStat();
+        QMutex mutex;
+        mutex.lock();
+        renderer->mRenderingFinished.wait(&mutex);
+        mutex.unlock();
+        mStats <<
+          "#AvgPvsRenderError\n" <<renderer->mPvsStat.GetAvgError()<<endl<<
+          "#MaxPvsRenderError\n" <<renderer->mPvsStat.GetMaxError()<<endl<<
+          "#ErrorFreeFrames\n" <<renderer->mPvsStat.GetErrorFreeFrames()<<endl;
+  }
+}
+
 bool
 RssPreprocessor::ComputeVisibility()
 {
+  
+  connect(this, SIGNAL(EvalPvsStat()), renderer, SLOT(EvalPvsStat()) );
+  
   cout<<"Rss Preprocessor started\n"<<flush;
   cout<<"Memory/ray "<<sizeof(VssRay)+sizeof(RssTreeNode::RayInfo)<<endl;
-  mSceneGraph->CollectObjects(&mObjects);
+
+  Randomize(0);
+  
         
   long startTime = GetTime();
@@ -443 +521 @@
 
   RssTree *rssTree = NULL;
-
-  if (mLoadInitialSamples)
-  {
-          cout << "Loading samples from file ... ";
-          LoadSamples(mVssRays, mObjects);
-          cout << "finished\n" << endl;
-  }
-  else
-  {
+  int rssPass = 0;
+  int rssSamples = 0;
+
+#if 0
+  if (mLoadInitialSamples) {
+        cout << "Loading samples from file ... ";
+        LoadSamples(mVssRays, mObjects);
+        cout << "finished\n" << endl;
+  } else {
+#endif
         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);
+          int s = mInitialSamples;
+          
+          for (int k=0; k < s; k++) {
+                
+                
                 //Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
                 Vector3 viewpoint; 
-                mViewCellsManager->GetViewPoint(viewpoint);
-                Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
-                                
+                //              mViewCellsManager->GetViewPoint(viewpoint);
+                viewpoint = GetViewpoint(mViewSpaceBox);
+                Vector3 direction = InitialGetDirection(viewpoint, mViewSpaceBox);
+                
                 sampleContributions = CastRay(viewpoint, direction, mVssRays);
-                                
-                                
+                
+                
                 //-- CORR matt: put block inside loop
                 if (sampleContributions) {
-                        passContributingSamples ++;
-                        passSampleContributions += sampleContributions;
+                  passContributingSamples ++;
+                  passSampleContributions += sampleContributions;
                 }
                 passSamples++;
                 totalSamples++;
-                }
-            
-                        
-                float avgRayContrib = (passContributingSamples > 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<<"%)" 
-                        << "avgcontrib=" << avgRayContrib << endl;
-                        
-                mStats <<
+          
+          cout << "#Pass " << mPass << " : t = " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+          cout << "#TotalSamples=" << totalSamples/1000 
+                   << "k   #SampleContributions=" << passSampleContributions << " (" 
+                   << 100*passContributingSamples/(float)passSamples<<"%)" 
+                   << "avgcontrib=" << avgRayContrib << endl;
+          
+          mPass++;
+        }
+          
+        cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
+        cout << "#totalRayStackSize=" << (int)mVssRays.size() << endl <<flush;
+        
+        rssSamples += mVssRays.size();
+        
+        
+        if (mExportRays) {
+          char filename[64];
+          sprintf(filename, "rss-rays-initial.x3d");
+          ExportRays(filename, mVssRays, mExportNumRays);
+        }
+        
+        if (mStoreInitialSamples)
+          {
+                cout << "Writing samples to file ... ";
+                WriteSamples(mVssRays);
+                cout << "finished\n" << endl;
+          }
+        
+        if (mUseViewcells) {
+          // construct view cells
+          mViewCellsManager->Construct(mObjects, mVssRays);
+          // evaluate contributions of the intitial rays
+          mViewCellsManager->ComputeSampleContributions(mVssRays);
+          
+          
+          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 <<
-                "#AvgRayContrib\n" << avgRayContrib << endl;
-                
-                mPass++;
-
-        }
-  
-        cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
-  }
-
-  cout << "#totalRayStackSize=" << (int)mVssRays.size() << endl <<flush;
-  
-
-  if (mExportRays) {
-        char filename[64];
-        sprintf(filename, "rss-rays-initial.x3d");
-        ExportRays(filename, mVssRays, mExportNumRays);
-  }
-  
-  if (mStoreInitialSamples)
-  {
-          cout << "Writing samples to file ... ";
-          WriteSamples(mVssRays);
-          cout << "finished\n" << endl;
-  }
-
-  if (mUseViewcells) {
-        // construct view cells
-        mViewCellsManager->Construct(mObjects, mVssRays);
-        
+                "#RssPass\n" <<rssPass<<endl<<
+                "#Time\n" << TimeDiff(startTime, GetTime())*1e-3<<endl<<
+                "#TotalSamples\n" <<totalSamples<<endl<<
+                "#RssSamples\n" <<rssSamples<<endl;
+          
+        ComputeRenderError();
+
+        mVssRays.PrintStatistics(mStats);
+        mViewCellsManager->PrintPvsStatistics(mStats);
+
         VssRayContainer selectedRays;
         int desired = Max(
                                           mViewCellsManager->GetPostProcessSamples(),
                                           mViewCellsManager->GetVisualizationSamples());
-        float p = desired/(float)mVssRays.size();
-        
-        for (int i=0; i < mVssRays.size(); i++) {
-          if (Random(1.0f) < p)
-                selectedRays.push_back(mVssRays[i]);
-        }
+
+        mVssRays.SelectRays(desired, selectedRays);
+        
         //-- post process view cells
         mViewCellsManager->PostProcess(mObjects, selectedRays);
@@ -547 +633 @@
 
 
+  
+  rssPass++;
+  
   rssTree = new RssTree;
-        
+  rssTree->SetPass(mPass);
+  
   if (mUseImportanceSampling) {
         
@@ -555 +645 @@
         else
           rssTree->Construct(mVssRays, NULL);
+
+        rssTree->stat.Print(mStats);
         cout<<"RssTree root PVS size = "<<rssTree->GetRootPvsSize()<<endl;
         
@@ -572 +664 @@
   // viewcells->UpdatePVS(newVssRays);
   
-  int rssPass = 0;
-  int rssSamples = 0;
+
   while (1) {
         int num = mRssSamplesPerPass;
         SimpleRayContainer rays;
         VssRayContainer vssRays;
-                
+
+
         if (!mUseImportanceSampling) {
           for (int j=0; j < num; j++) {
@@ -597 +689 @@
 
         
-        totalSamples+=num;
-        rssSamples+=num;
-
+        totalSamples += rays.size();
+        rssSamples += vssRays.size();
+
+        
+        
         mStats <<
           "#Pass\n" <<mPass<<endl<<
@@ -607 +701 @@
           "#RssSamples\n" <<rssSamples<<endl;
 
+
         if (mUseViewcells) {
           
@@ -616 +711 @@
         }
 
-        // add rays to the tree after the viewcells have been cast to have their contributions
-        // already when adding into the tree
-        // do not add those rays which have too low or no contribution....
-        
-        if (mUseImportanceSampling) {
-          rssTree->AddRays(vssRays);
-          
-          if (0) {
-                cout<<"############# Rss PVS STAT ##################\n";
-                cout<<"#AVG_RSS_PVS\n"<<rssTree->GetAvgPvsSize()<<endl;
-                cout<<"#RSS_ROOT_PVS\n"<<rssTree->GetRootPvsSize()<<endl;
-          }
-          
-          if (mUpdateSubdivision) {
-                int subdivided = rssTree->UpdateSubdivision();
-                cout<<"subdivided leafs = "<<subdivided<<endl;
-                cout<<"#total leaves = "<<rssTree->stat.Leaves()<<endl;
-          }
-        }
-        
+
+        ComputeRenderError();
+        
+        // epxort rays before adding them to the tree -> some of them can be deleted
+
         if (mExportRays) {
           char filename[64];
@@ -644 +724 @@
           
           ExportRays(filename, vssRays, mExportNumRays);
-        }
+
+          // now export all contributing rays
+          VssRayContainer contributingRays;
+          vssRays.GetContributingRays(contributingRays, mPass);
+          mStats<<"#NUM_CONTRIBUTING_RAYS\n"<<contributingRays.size()<<endl;
+          sprintf(filename, "rss-crays-%04d.x3d", rssPass);
+          ExportRays(filename, contributingRays, mExportNumRays);
+        }
+
+        
+        // add rays to the tree after the viewcells have been cast to have their contributions
+        // already when adding into the tree
+        // do not add those rays which have too low or no contribution....
+        
+        if (mUseImportanceSampling) {
+          rssTree->AddRays(vssRays);
+          
+          if (0) {
+                cout<<"############# Rss PVS STAT ##################\n";
+                cout<<"#AVG_RSS_PVS\n"<<rssTree->GetAvgPvsSize()<<endl;
+                cout<<"#RSS_ROOT_PVS\n"<<rssTree->GetRootPvsSize()<<endl;
+          }
+          
+          
+          if (mUpdateSubdivision) {
+                int updatePasses = 1;
+                if (mPass % updatePasses == 0) {
+                  int subdivided = rssTree->UpdateSubdivision();
+                  cout<<"subdivided leafs = "<<subdivided<<endl;
+                  cout<<"#total leaves = "<<rssTree->stat.Leaves()<<endl;
+                }
+          }
+        }
+        
 
         
@@ -653 +766 @@
         }
 
+        
+        if (!mUseImportanceSampling) 
+          CLEAR_CONTAINER(vssRays);
+
         if (totalSamples >= mRssSamples + mInitialSamples)
           break;
@@ -659 +776 @@
         rssPass++;
         mPass++;
+        rssTree->SetPass(mPass);
   }
   

trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.h

@@ -4 +4 @@
 #include <fstream>
 using namespace std;
+
 
 #include "Preprocessor.h"
@@ -14 +15 @@
     sampling of view space */
 class RssPreprocessor : public Preprocessor {
+  
 public:
   int mPass;
@@ -29 +31 @@
   bool mUseViewcells;
   bool mUpdateSubdivision;
+
+  bool mObjectBasedSampling;
+  bool mDirectionalSampling;
   
   AxisAlignedBox3 *mViewSpaceBox;
@@ -34 +39 @@
   ofstream mStats;
 
-  ObjectContainer mObjects;
 
   // rays cast during the processing
@@ -56 +60 @@
                            AxisAlignedBox3 *viewSpaceBox
                            );
-        
+
+  Vector3
+  InitialGetDirection(const Vector3 &viewpoint,
+                                          AxisAlignedBox3 *viewSpaceBox
+                                          );
+
   void
   SetupRay(Ray &ray, 
@@ -118 +127 @@
                         );
 
+
+  void
+  ComputeRenderError();
+
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.cpp

@@ -1 @@
-
 // ================================================================
 // $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
@@ -25 +24 @@
 #include "Ray.h"
 #include "Containers.h"
+#include "ViewCell.h"
+#include "Exporter.h"
+#include "Preprocessor.h"
+#include "SceneGraph.h"
 
 #define DEBUG_SPLIT_COST 0
@@ -64 +67 @@
 }
 
+inline void
+AddViewcells2Pvs(const ViewCellContainer &viewcells,
+                                 const int side,
+                                 int &viewcellsBack,
+                                 int &viewcellsFront)
+{
+  ViewCellContainer::const_iterator it = viewcells.begin();
+  
+  for (; it != viewcells.end(); ++it) {
+        ViewCell *viewcell = *it;
+        if (side <= 0) {
+          if (!viewcell->Mailed() && !viewcell->Mailed(2)) {
+                viewcellsBack++;
+                if (viewcell->Mailed(1))
+                  viewcell->Mail(2);
+                else
+                  viewcell->Mail();
+          }
+        }
+        
+        if (side >= 0) {
+          if (!viewcell->Mailed(1) && !viewcell->Mailed(2)) {
+                viewcellsFront++;
+                if (viewcell->Mailed())
+                  viewcell->Mail(2);
+                else
+                  viewcell->Mail(1);
+          }
+        }
+  }
+}
+
+
 // Constructor
 RssTree::RssTree()
@@ -85 +121 @@
   environment->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
   
-  
-
-
-  float refDirAngle;
-  environment->GetFloatValue("RssTree.refDirAngle", refDirAngle);
+  environment->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
+
+
   
   environment->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
@@ -125 +159 @@
   environment->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);
 
+  environment->GetBoolValue("RssTree.importanceBasedCost", mImportanceBasedCost);
+  
+  environment->GetIntValue("RssTree.maxRays", mMaxRays);
+
   root = NULL;
   
@@ -143 +181 @@
 RssStatistics::Print(ostream &app) const
 {
-  app << "===== RssTree statistics ===============\n";
+  app << "###### RssTree statistics ######\n";
 
   app << "#N_RAYS ( Number of rays )\n"
@@ -204 +242 @@
       << Time() << " \n";
 
-  app << "===== END OF RssTree statistics ==========\n";
+  app << "###### END OF RssTree statistics ######\n";
 
 }
@@ -210 +248 @@
 
 void
-RssTreeLeaf::UpdatePvsSize()
-{
-  if (!mValidPvs) {
+RssTree::UpdatePvsSize(RssTreeLeaf *leaf)
+{
+  if (!leaf->mValidPvs) {
         Intersectable::NewMail();
         int pvsSize = 0;
-        for(RssTreeNode::RayInfoContainer::iterator ri = rays.begin();
-                ri != rays.end();
+        for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                ri != leaf->rays.end();
                 ri++)
           if ((*ri).mRay->IsActive()) {
                 Intersectable *object;
-#if BIDIRECTIONAL_RAY
-                object = (*ri).mRay->mOriginObject;
+                object = (*ri).GetObject();
                 if (object && !object->Mailed()) {
                   pvsSize++;
                   object->Mail();
                 }
-#endif
-                object = (*ri).mRay->mTerminationObject;
-                if (object && !object->Mailed()) {
-                  pvsSize++;
-                  object->Mail();
-                }
-          }
-        mPvsSize = pvsSize;
-        mValidPvs = true;
-
-        ComputeEntropyImportance();
+          }
+        leaf->SetPvsSize(pvsSize);
+        
+        ComputeImportance(leaf);
   }
 }
@@ -248 +278 @@
   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)
+  cerr<<"Clip not reimplented yet...Quiting\n";
+  exit(1);
+  ray.Init(rayInfo.GetOrigin(), rayInfo.GetDir(), Ray::LINE_SEGMENT);
+
+  if (!box.ComputeMinMaxT(ray, &tmin, &tmax) || tmin>=tmax)
         return false;
   
@@ -281 +313 @@
   if (root)
     delete root;
-        
+
+  
   root = new RssTreeLeaf(NULL, rays.size());
   // first construct a leaf that will get subdivide
@@ -315 +348 @@
                                         );
   }
-        
-        
+  
+  // make the z axis (unused) a unit size
+  // important for volume computation
+  dirBBox.SetMin(2, 0.0f);
+  dirBBox.SetMax(2, 1.0f);
+  
   if ( forcedBoundingBox ) 
         bbox = *forcedBoundingBox;
@@ -324 +361 @@
 
   stat.rays = leaf->rays.size();
-  leaf->UpdatePvsSize();
-  leaf->ComputeEntropyImportance();
+  UpdatePvsSize(leaf);
 
   stat.initialPvsSize = leaf->GetPvsSize();
@@ -338 +374 @@
   
   stat.Stop();
-
   stat.Print(cout);
   cout<<"#Total memory="<<GetMemUsage()<<endl;
 
+  // this rotine also updates importances etc...
 }
 
@@ -387 +423 @@
           if (!node->IsLeaf()) {
                 subdivided++;
+
+                
                 RssTreeInterior *interior = (RssTreeInterior *) node;
                 // push the children on the stack
@@ -482 +520 @@
                                          RssTreeLeaf *leaf,
                                          const AxisAlignedBox3 &box,
-                                         float &position,
-                                         int &raysBack,
-                                         int &raysFront,
-                                         int &pvsBack,
-                                         int &pvsFront
+                                         SplitInfo &info
                                          )
 {
 
-  int minDirDepth = 6;
-  int axis;
-  float costRatio;
-        
-  costRatio = BestCostRatio(leaf,
-                                                        axis,
-                                                        position,
-                                                        raysBack,
-                                                        raysFront,
-                                                        pvsBack,
-                                                        pvsFront
-                                                        );
+  BestCostRatio(leaf,
+                                info);
+  
 #if DEBUG_SPLIT_COST
-  cout<<axis<<" r="<<costRatio<<endl;
+  Debug<<"Split Info:"<<endl;
+  Debug<<"axis="<<info.axis<<" ratio="<<info.costRatio<<endl;
+  Debug<<"viewcells="<<info.viewCells<<
+        " viewcells back="<<info.viewCellsBack<<
+        " viewcells back="<<info.viewCellsFront<<endl;
 #endif
-  if (costRatio > termMaxCostRatio) {
+  
+  if (info.costRatio > termMaxCostRatio) {
         //              cout<<"Too big cost ratio "<<costRatio<<endl;
         stat.maxCostRatioNodes++;
@@ -516 +546 @@
         " rays="<<leaf->rays.size()<<
         " rc="<<leaf->GetAvgRayContribution()<<
-        " axis="<<axis<<endl;
+        " axis="<<info.axis<<endl;
 #endif
-        
-  return axis;
-}
-
-
-float
+  
+  return info.axis;
+}
+
+
+void
 RssTree::GetCostRatio(
                                           RssTreeLeaf *leaf,
-                                          const int axis,
-                                          const float position,
-                                          const int raysBack,
-                                          const int raysFront,
-                                          const int pvsBack,
-                                          const int pvsFront
+                                          SplitInfo &info
                                           )
 {
-  bool costImportance = true;
-
-  float ratio;
+
   AxisAlignedBox3 box;
   float minBox, maxBox;
 
-  if (axis < 3) {
+  if (info.axis < 3) {
         box = GetBBox(leaf);
-        minBox = box.Min(axis);
-        maxBox = box.Max(axis);
+        minBox = box.Min(info.axis);
+        maxBox = box.Max(info.axis);
   }     else {
         box = GetDirBBox(leaf);
-        minBox = box.Min(axis-3);
-        maxBox = box.Max(axis-3);
+        minBox = box.Min(info.axis-3);
+        maxBox = box.Max(info.axis-3);
   }
         
@@ -554 +577 @@
   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;
+  if (!mImportanceBasedCost) {
+        const int costMethod = 0;
+        
+        switch (costMethod) {
+        case 0: {
+          //            float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+          float sum = info.pvsBack*(info.position - minBox) + info.pvsFront*(maxBox - info.position);
+          float newCost = ct_div_ci + sum/sizeBox;
+          float oldCost = pvsSize;
+          info.costRatio = newCost/oldCost;
+          break;
+        }
+        case 1: {
+          float newContrib =
+                info.contributionBack/(info.position - minBox) + 
+                +
+                info.contributionFront/(maxBox - info.position);
+          float oldContrib = info.contribution/sizeBox;
+          info.costRatio = oldContrib/newContrib;
+          break;
+        } 
+        case 2: {
+          float sum =
+                info.viewCellsBack*(info.position - minBox) +
+                info.viewCellsFront*(maxBox - info.position);
+          float newCost = ct_div_ci + sum/sizeBox;
+          float oldCost = info.viewCells;
+          info.costRatio = newCost/oldCost;
+          break;
+        }
+        case 3: {
+          float newCost = info.raysBack*info.pvsBack  + info.raysFront*info.pvsFront;
+          float oldCost = leaf->rays.size()*pvsSize;
+          info.costRatio = newCost/oldCost;
+        }
+        }
   } else {
+        const int costMethod = 1;
         // importance based cost
-#if 0
-        float newContrib =
-          ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
-           (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
-                
-        //                      float newContrib =
-        //                              sqr(pvsBack/(raysBack + Limits::Small)) +
-        //                              sqr(pvsFront/(raysFront + Limits::Small));
-        float oldContrib = sqr(leaf->GetAvgRayContribution());
-        ratio = oldContrib/newContrib;
-#else
-#if 1
-        float newCost = raysBack*pvsBack  + raysFront*pvsFront;
-        float oldCost = leaf->rays.size()*pvsSize;
-        ratio = newCost/oldCost;
-#else
-#if 0
-        float newCost = (pvsBack  + pvsFront)*0.5f;
-        float oldCost = pvsSize;
-        ratio = newCost/oldCost;
-#else
-        float newCost = abs(raysBack  - raysFront);
-        float oldCost = leaf->rays.size();
-        ratio = newCost/oldCost;
-#endif
-#endif
-#endif
-  }
-        
-  return ratio;
+        switch (costMethod) {
+        case 0: {
+          break;
+        }
+        case 1: {
+          float newContrib =
+                sqr(info.pvsBack/(info.raysBack + Limits::Small)) +
+                sqr(info.pvsFront/(info.raysFront + Limits::Small));
+          float oldContrib = sqr(leaf->GetAvgRayContribution());
+          info.costRatio = oldContrib/newContrib;
+          break;
+        }
+        case 2: {
+          float newCost = (info.pvsBack  + info.pvsFront)*0.5f;
+          float oldCost = pvsSize;
+          info.costRatio = newCost/oldCost;
+          break;
+        }
+        case 3: {
+          float newCost = abs(info.raysBack  - info.raysFront);
+          float oldCost = leaf->rays.size();
+          info.costRatio = newCost/oldCost;
+          break;
+        }
+  }
+  }
 }
                                                         
 
-float
+void
 RssTree::EvalCostRatio(
                                            RssTreeLeaf *leaf,
-                                           const int axis,
-                                           const float position,
-                                           int &raysBack,
-                                           int &raysFront,
-                                           int &pvsBack,
-                                           int &pvsFront
+                                           SplitInfo &info
                                            )
 {
-  raysBack = 0;
-  raysFront = 0;
-  pvsFront = 0;
-  pvsBack = 0;
-
-        
+  info.raysBack = 0;
+  info.raysFront = 0;
+  info.pvsFront = 0;
+  info.pvsBack = 0;
+  info.viewCells = 0;
+  info.viewCellsFront = 0;
+  info.viewCellsBack = 0;
+
+
+  float sumWeights = Limits::Small;
+  float sumWeightsBack = Limits::Small;
+  float sumWeightsFront = Limits::Small;
+  
+  float sumContribution = 0.0f;
+  float sumContributionBack = 0.0f;
+  float sumContributionFront = 0.0f;
+
+  Intersectable::NewMail();
+
+  // count the numebr of viewcells first
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+          ri != leaf->rays.end();
+          ri++)
+        if ((*ri).mRay->IsActive()) {
+          ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
+          for (; it != (*ri).mRay->mViewCells.end(); ++it) {
+                if (!(*it)->Mailed()) {
+                  (*it)->Mail();
+                  info.viewCells++;
+                }
+          }
+        }
+  
   Intersectable::NewMail(3);
-        
-  if (axis <= RssTreeNode::SPLIT_Z) {
-    // this is the main ray classification loop!
-    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                ri != leaf->rays.end();
-                ri++)
-      if ((*ri).mRay->IsActive()) {
-                                
-                // determine the side of this ray with respect to the plane
-                int side = (*ri).ComputeRaySide(axis, position);
-                //                              (*ri).mRay->mSide = side;
-                                
-                if (side <= 0)
-                  raysBack++;
-                                
-                if (side >= 0)
-                  raysFront++;
-                                
-                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
-      }
-                
-  } else {
-                
-        // directional split
-    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                ri != leaf->rays.end();
-                ri++)
-      if ((*ri).mRay->IsActive()) {
-                                
-                // determine the side of this ray with respect to the plane
-                int side;
-                if ((*ri).mRay->GetDirParametrization(axis - 3) <= position)
-                  side = -1;
-                else
-                  side = 1;
-                
-                if (side <= 0)
-                  raysBack++;
-                                
-                if (side >= 0)
-                  raysFront++;
-
-                //                              (*ri).mRay->mSide = side;
-                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
-                                
-      }
-  }
-
-  float ratio = GetCostRatio(
-                                                         leaf,
-                                                         axis,
-                                                         position,
-                                                         raysBack,
-                                                         raysFront,
-                                                         pvsBack,
-                                                         pvsFront);
-
+  
+  // this is the main ray classification loop!
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+          ri != leaf->rays.end();
+          ri++)
+        if ((*ri).mRay->IsActive()) {
+          int side;
+          
+          
+          // determine the side of this ray with respect to the plane
+          side = (*ri).ComputeRaySide(info.axis, info.position);
+          
+          float weight, contribution;
+          
+          GetRayContribution(*ri, weight, contribution);
+          sumWeights += weight;
+          sumContribution += contribution;
+          
+          if (side <= 0) {
+                info.raysBack++;
+                sumWeightsBack += weight;
+                sumContributionBack += contribution;
+          }
+          
+          if (side >= 0) {
+                info.raysFront++;
+                sumWeightsFront += weight;
+                sumContributionFront += contribution;
+          }
+          
+          AddObject2Pvs((*ri).GetObject(), side, info.pvsBack, info.pvsFront);
+          AddViewcells2Pvs((*ri).mRay->mViewCells,
+                                           side,
+                                           info.viewCellsBack,
+                                           info.viewCellsFront);
+        }
+
+  info.contribution = sumContribution/sumWeights;
+  info.contributionBack = sumContributionBack/sumWeightsBack;
+  info.contributionFront = sumContributionFront/sumWeightsFront;
+  
+  GetCostRatio(
+                           leaf,
+                           info);
+  
   //    cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
   //  float oldCost = leaf->rays.size();
-
+  
   //    cout<<"ratio="<<ratio<<endl;
-
-  return ratio;
-}
-
-float
+}
+
+void
 RssTree::BestCostRatio(
                                            RssTreeLeaf *leaf,
-                                           int &axis,
-                                           float &position,
-                                           int &raysBack,
-                                           int &raysFront,
-                                           int &pvsBack,
-                                           int &pvsFront
+                                           SplitInfo &info
                                            )
 {
-  int nRaysBack[6], nRaysFront[6];
-  int nPvsBack[6], nPvsFront[6];
-  float nPosition[6];
-  float nCostRatio[6];
+  SplitInfo nInfo[5];
   int bestAxis = -1;
-        
+  
   AxisAlignedBox3 sBox = GetBBox(leaf);
   AxisAlignedBox3 dBox = GetDirBBox(leaf);
@@ -700 +750 @@
   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++)
+  for (int axis = 0; axis < 5; axis++)
         if (
                 (axis < 3 && (leaf->depth < mDirSplitDepth ||  mInterleaveDirSplits)) ||
                 (axis >= 3 && (leaf->depth >= mDirSplitDepth))
                 ) {
+          nInfo[axis].axis = axis;
           if (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis) {
                 
                 if (splitType == ESplitRegular) {
                   if (axis < 3)
-                        nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
+                        nInfo[axis].position = (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]
-                                                                                   );
+                        nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
+                  EvalCostRatio(leaf,
+                                                nInfo[axis]);
                 } else
                   if (splitType == ESplitHeuristic) {
-                        nCostRatio[axis] = EvalCostRatioHeuristic(
-                                                                                                          leaf,
-                                                                                                          axis,
-                                                                                                          nPosition[axis],
-                                                                                                          nRaysBack[axis],
-                                                                                                          nRaysFront[axis],
-                                                                                                          nPvsBack[axis],
-                                                                                                          nPvsFront[axis]);
+                        EvalCostRatioHeuristic(
+                                                                   leaf,
+                                                                   nInfo[axis]
+                                                                   );
                   } else
                         if (splitType == ESplitHybrid) {
                           if (leaf->depth > 7)
-                                nCostRatio[axis] = EvalCostRatioHeuristic(
-                                                                                                                  leaf,
-                                                                                                                  axis,
-                                                                                                                  nPosition[axis],
-                                                                                                                  nRaysBack[axis],
-                                                                                                                  nRaysFront[axis],
-                                                                                                                  nPvsBack[axis],
-                                                                                                                  nPvsFront[axis]);
+                                EvalCostRatioHeuristic(
+                                                                           leaf,
+                                                                           nInfo[axis]
+                                                                           );
                           else {
                                 if (axis < 3)
-                                  nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
+                                  nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
                                 else
-                                  nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
+                                  nInfo[axis].position = (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]
-                                                                                                 );
+                                EvalCostRatio(leaf,
+                                                          nInfo[axis]
+                                                          );
                           }
                         } else {
@@ -768 +798 @@
                         }
                 
-                
                 if ( bestAxis == -1)
                   bestAxis = axis;
                 else
-                  if ( nCostRatio[axis] < nCostRatio[bestAxis] )
+                  if ( nInfo[axis].costRatio < nInfo[bestAxis].costRatio )
                         bestAxis = axis;
           }
         }
   
-  axis = bestAxis;
-  position = nPosition[bestAxis];
-
-  raysBack = nRaysBack[bestAxis];
-  raysFront = nRaysFront[bestAxis];
-
-  pvsBack = nPvsBack[bestAxis];
-  pvsFront = nPvsFront[bestAxis];
-        
-  return nCostRatio[bestAxis];
-}
-
-        
-float
+  info = nInfo[bestAxis];
+}
+
+        
+void
 RssTree::EvalCostRatioHeuristic(
                                                                 RssTreeLeaf *leaf,
-                                                                const int axis,
-                                                                float &bestPosition,
-                                                                int &raysBack,
-                                                                int &raysFront,
-                                                                int &pvsBack,
-                                                                int &pvsFront
+                                                                SplitInfo &info
                                                                 )
 {
@@ -804 +819 @@
   float minBox, maxBox;
         
-  if (axis < 3) {
+  if (info.axis < 3) {
         box = GetBBox(leaf);
-        minBox = box.Min(axis);
-        maxBox = box.Max(axis);
+        minBox = box.Min(info.axis);
+        maxBox = box.Max(info.axis);
   } else {
         box = GetDirBBox(leaf);
-        minBox = box.Min(axis-3);
-        maxBox = box.Max(axis-3);
-  }
-        
-  SortSplitCandidates(leaf, axis);
+        minBox = box.Min(info.axis-3);
+        maxBox = box.Max(info.axis-3);
+  }
+        
+  SortSplitCandidates(leaf, info.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();
+  
+  SplitInfo currInfo;
+  currInfo.axis = info.axis;
+
+  currInfo.raysBack = 0;
+  currInfo.raysFront = leaf->rays.size();
+  
+  currInfo.pvsBack = 0;
+  currInfo.pvsFront = leaf->GetPvsSize();
+
+ 
   float sizeBox = maxBox - minBox;
-        
+  
   float minBand = minBox + 0.1f*(maxBox - minBox);
   float maxBand = minBox + 0.9f*(maxBox - minBox);
         
-  float minRatio = 1e20f;
-        
+  // best cost ratio
+  info.costRatio = 1e20f;
+
+  currInfo.viewCells = 0;
+
   Intersectable::NewMail();
   // set all object as belonging to the fron pvs
@@ -835 +861 @@
           ri++)
         if ((*ri).mRay->IsActive()) {
-          Intersectable *object = (*ri).mRay->mTerminationObject;
+          Intersectable *object = (*ri).GetObject();
           if (object)
                 if (!object->Mailed()) {
@@ -842 +868 @@
                 } else
                   object->mCounter++;
-        }
-        
+          
+          // and do the same for all viewcells
+          ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
+          
+          for (; it != (*ri).mRay->mViewCells.end(); ++it) {
+                ViewCell *viewcell = *it;
+                if (!viewcell->Mailed()) {
+                  currInfo.viewCells++;
+                  viewcell->Mail();
+                  viewcell->mCounter = 1;
+                } else
+                  viewcell->mCounter++;
+          }
+        }
+
+  currInfo.viewCellsBack = 0;
+  currInfo.viewCellsFront = currInfo.viewCells;
+ 
   Intersectable::NewMail();
         
@@ -849 +891 @@
           ci < splitCandidates->end();
           ci++) {
-        VssRay *ray;
         switch ((*ci).type) {
         case SortableEntry::ERayMin: {
-          rr--;
-          rl++;
-          ray = (VssRay *) (*ci).data;
-          Intersectable *object = ray->mTerminationObject;
+          currInfo.raysFront--;
+          currInfo.raysBack++;
+          RssTreeNode::RayInfo *rayInfo = (RssTreeNode::RayInfo *) (*ci).data;
+          Intersectable *object = rayInfo->GetObject();
           if (object) {
                 if (!object->Mailed()) {
                   object->Mail();
-                  pl++;
+                  currInfo.pvsBack++;
                 }
                 if (--object->mCounter == 0)
-                  pr--;
+                  currInfo.pvsFront--;
+          }
+          ViewCellContainer::const_iterator it = rayInfo->mRay->mViewCells.begin();
+          for (; it != rayInfo->mRay->mViewCells.end(); ++it) {
+                ViewCell *viewcell = *it;
+                if (!viewcell->Mailed()) {
+                  viewcell->Mail();
+                  currInfo.viewCellsBack++;
+                }
+                if (--viewcell->mCounter == 0)
+                  currInfo.viewCellsFront--;
           }
           break;
@@ -871 +922 @@
         
         if (position > minBand && position < maxBand) {
-                        
-          float ratio = GetCostRatio(
-                                                                 leaf,
-                                                                 axis,
-                                                                 position,
-                                                                 rl,
-                                                                 rr,
-                                                                 pl,
-                                                                 pr);
-                        
+          currInfo.position = position;
+          
+          GetCostRatio(
+                                   leaf,
+                                   currInfo);
+          
                         
           //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
           //      cout<<"cost= "<<sum<<endl;
-                        
-          if (ratio < minRatio) {
-                minRatio = ratio;
-                bestPosition = position;
-                                
-                raysBack = rl;
-                raysFront = rr;
-                                
-                pvsBack = pl;
-                pvsFront = pr;
-                                
+          
+          if (currInfo.costRatio < info.costRatio) {
+                info = currInfo;
       }
     }
   }
-
+  
   
   //  cout<<"===================="<<endl;
   //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
   //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
-  return minRatio;
 }
 
@@ -933 +971 @@
           if (axis < 3) {
                 splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
-                                                                                                 (*ri).ExtrapOrigin(axis),
-                                                                                                 (void *)(*ri).mRay)
+                                                                                                 (*ri).GetOrigin(axis),
+                                                                                                 (void *)&(*ri))
                                                                    );
           } else {
-                float pos = (*ri).mRay->GetDirParametrization(axis-3);
+                float pos = (*ri).GetDirParametrization(axis-3);
                 splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
                                                                                                  pos,
-                                                                                                 (void *)(*ri).mRay)
+                                                                                                 (void *)&(*ri))
                                                                    );
           }
@@ -1012 +1050 @@
   }
         
-  float position;
-        
-  // first count ray sides
-  int raysBack;
-  int raysFront;
-  int pvsBack;
-  int pvsFront;
+  SplitInfo info;
         
   // select subdivision axis
-  int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
-  Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
+  int axis = SelectPlane( leaf,
+                                                  box,
+                                                  info
+                                                  );
+  //  Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
   
   if (axis == -1) {
@@ -1035 +1070 @@
 
   node->axis = axis;
-  node->position = position;
+  node->position = info.position;
   node->bbox = box;
   node->dirBBox = GetDirBBox(leaf);
@@ -1041 +1076 @@
   backBBox = box;
   frontBBox = box;
-
-  RssTreeLeaf *back = new RssTreeLeaf(node, raysBack);
-  RssTreeLeaf *front = new RssTreeLeaf(node, raysFront);
-
+  
+  RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
+  RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);
+
+  // update halton generator
+  back->halton.index = leaf->halton.index;
+  front->halton.index = leaf->halton.index;
+  
   // replace a link from node's parent
   if (  leaf->parent )
@@ -1052 +1091 @@
         
   if (axis <= RssTreeNode::SPLIT_Z) {
-        backBBox.SetMax(axis, position);
-    frontBBox.SetMin(axis, position);
-                
-        for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                ri != leaf->rays.end();
-                ri++) {
-      if ((*ri).mRay->IsActive()) {
-                                
-                // first unref ray from the former leaf
-                (*ri).mRay->Unref();
-
-                // Debug << "computed t: " << (*ri).mRay->mT << endl;
-                // determine the side of this ray with respect to the plane
-                int side = node->ComputeRaySide(*ri);
-                
-                if (side == 0) {
-                  if ((*ri).mRay->HasPosDir(axis)) {
-                        back->AddRay(*ri);
-                        front->AddRay(*ri);
-                  } else {
-                        back->AddRay(*ri);
-                        front->AddRay(*ri);
-                  }
-                } else
-                  if (side == 1) 
-                        front->AddRay(*ri);
-                  else
-                        back->AddRay(*ri);
-      } else
-                (*ri).mRay->Unref();
-    }
-  } else {
-    // rays front/back
-    
-    
-    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                ri != leaf->rays.end();
-                ri++) {
-      if ((*ri).mRay->IsActive()) {
-                // first unref ray from the former leaf
-                (*ri).mRay->Unref();
-
-                int side;
-                if ((*ri).mRay->GetDirParametrization(axis - 3) <= position)
-                  side = -1;
-                else
-                  side = 1;
-                                
-                if (side == 1)
-                  front->AddRay(*ri);
-                else
-                  back->AddRay(*ri);
-                                
-      } else
-                (*ri).mRay->Unref();
-    }
-  }
-
+        backBBox.SetMax(axis, info.position);
+    frontBBox.SetMin(axis, info.position);
+  }
+
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+          ri != leaf->rays.end();
+          ri++) {
+        if ((*ri).mRay->IsActive()) {
+          
+          // first unref ray from the former leaf
+          (*ri).mRay->Unref();
+          
+          // Debug << "computed t: " << (*ri).mRay->mT << endl;
+          // determine the side of this ray with respect to the plane
+          int side = node->ComputeRaySide(*ri);
+          
+          if (side == 1) 
+                front->AddRay(*ri);
+          else
+                back->AddRay(*ri);
+        } else
+          (*ri).mRay->Unref();
+  }
+
+  // distribute the total number of rays according to the distribution
+  // of rays which remained
+  
+  
+  //  front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
+  //  back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
+  
 #if 0
   front->SetPvsSize(pvsFront);
@@ -1117 +1129 @@
   back->ComputeEntropyImportance();
 #else
-  front->UpdatePvsSize();
-  back->UpdatePvsSize();
+  UpdatePvsSize(front);
+  UpdatePvsSize(back);
 #endif
 
@@ -1124 +1136 @@
   // update stats
   stat.rayRefs -= (int)leaf->rays.size();
-  stat.rayRefs += raysBack + raysFront;
+  stat.rayRefs += info.raysBack + info.raysFront;
 
   
@@ -1221 +1233 @@
 
 void
-RssTree::UpdateRays(VssRayContainer &remove,
+RssTree::UpdateRays(
+                                        VssRayContainer &remove,
                                         VssRayContainer &add
                                         )
@@ -1257 +1270 @@
           AddRay(info);
   }
-}
+
+  stat.rayRefs += add.size() - remove.size();
+
+  UpdateTreeStatistics();
+  // check whether the tree should be prunned
+  if (stat.rayRefs > mMaxRays) {
+        PruneRays(mMaxRays);
+        //      UpdateTreeStatistics();
+  }
+  
+}
+
+  
 
 
@@ -1377 +1402 @@
   RssTreeInterior *in =  (RssTreeInterior *) data.node;
   
-  if (in->axis <= RssTreeNode::SPLIT_Z) {
-
-    // determine the side of this ray with respect to the plane
-    int side = in->ComputeRaySide(data.rayData
-                                                                  );
-    
-  
-    if (side == 0) {
-      if (data.rayData.mRay->HasPosDir(in->axis)) {
-                tstack.push(RayTraversalData(in->back,
-                                                                         data.rayData)
-                                        );
-                                
-                tstack.push(RayTraversalData(in->front,
-                                                                         data.rayData)
-                                        );
-        
-      } else {
-                tstack.push(RayTraversalData(in->back,
-                                                                         data.rayData
-                                                                         )
-                                        );
-                                
-                tstack.push(RayTraversalData(in->front,
-                                                                         data.rayData)
-                                        );
-                                
-                                
-      }
-    } else
-      if (side == 1)
-                tstack.push(RayTraversalData(in->front, data.rayData));
-      else
-                tstack.push(RayTraversalData(in->back, data.rayData));
-  } 
-  else {
-    // directional split
-        if (data.rayData.mRay->GetDirParametrization(in->axis - 3) <= in->position)
-          tstack.push(RayTraversalData(in->back, data.rayData));
-        else
-          tstack.push(RayTraversalData(in->front, data.rayData));
-  }
+  // determine the side of this ray with respect to the plane
+  int side = in->ComputeRaySide(data.rayData
+                                                                );
+  
+  if (side == 1)
+        tstack.push(RayTraversalData(in->front, data.rayData));
+  else
+        tstack.push(RayTraversalData(in->back, data.rayData));
+  
 }
 
@@ -1559 +1551 @@
                 if ((*ri).mRay->IsActive()) {
                   Intersectable *object;
-#if BIDIRECTIONAL_RAY
-                  object = (*ri).mRay->mOriginObject;
-                  if (object && !object->Mailed()) {
-                        pvsSize++;
-                        object->Mail();
-                  }
-#endif
-                  object = (*ri).mRay->mTerminationObject;
+                  object = (*ri).GetObject();
                   if (object && !object->Mailed()) {
                         pvsSize++;
@@ -1612 +1597 @@
                 if ((*ri).mRay->IsActive()) {
                   Intersectable *object;
-                  object = (*ri).mRay->mTerminationObject;
+                  object = (*ri).GetObject();
                   if (object && !object->Mailed()) {
                         pvs.push_back(object);
@@ -1637 +1622 @@
 
 void
-RssTree::GetTreeStatistics(
-                                                   float &avgPvsSize,
-                                                   float &avgRays,
-                                                   float &avgRayContribution,
-                                                   float &avgPvsEntropy,
-                                                   float &avgRayLengthEntropy,
-                                                   float &avgImportance
-                                                   )
+RssTree::UpdateTreeStatistics()
 {
   stack<RssTreeNode *> tstack;
@@ -1651 +1629 @@
   float sumPvsSize = 0.0f;
   float sumRayContribution = 0.0f;
+  float sumWeightedRayContribution = 0.0f;
   float sumImportance = 0.0f;
   float sumPvsEntropy = 0.0f;
@@ -1665 +1644 @@
           leaves++;
           RssTreeLeaf *leaf = (RssTreeLeaf *)node;
-          leaf->UpdatePvsSize();
+          UpdatePvsSize(leaf);
           
           sumPvsSize += leaf->GetPvsSize();
           sumRayContribution += leaf->GetAvgRayContribution();
-          sumPvsEntropy += leaf->mPvsEntropy;
-          sumRayLengthEntropy += leaf->mRayLengthEntropy;
+          
+
+          RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
+          for (;it != leaf->rays.end(); ++it) {
+                float weight, contribution;
+                GetRayContribution(*it, weight, contribution);
+                sumWeightedRayContribution += weight*contribution;
+          }
+          
+          //      sumPvsEntropy += leaf->mPvsEntropy;
+          //      sumRayLengthEntropy += leaf->mRayLengthEntropy;
           sumRays += leaf->rays.size();
           
           float imp = leaf->GetImportance();
           
-          if (imp > 1.0f)
-                cout<<"warning imp > 1.0f:"<<imp<<endl;
+          //      if (imp > 1.0f)
+          //            cout<<"warning imp > 1.0f:"<<imp<<endl;
           
           sumImportance += imp;
@@ -1689 +1677 @@
   }
   
-  avgPvsSize = sumPvsSize/(float)leaves;
-  avgRays = sumRays/(float)leaves;
-  avgRayContribution = sumRayContribution/(float)leaves;
-  avgPvsEntropy = sumPvsEntropy/(float)leaves;
-  avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
-  avgImportance = sumImportance/(float)leaves;
-}
+  stat.avgPvsSize = sumPvsSize/(float)leaves;
+  stat.avgRays = sumRays/(float)leaves;
+  stat.avgRayContribution = sumRayContribution/(float)leaves;
+  //  avgPvsEntropy = sumPvsEntropy/(float)leaves;
+  //  avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
+  stat.avgImportance = sumImportance/(float)leaves;
+  stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
+  stat.rayRefs = (int)sumRays;
+}
+
 
 
 int
-RssTree::GenerateRays(const float ratioPerLeaf,
+RssTree::GenerateRays(
+                                          const float ratioPerLeaf,
                                           SimpleRayContainer &rays)
 {
@@ -1777 +1769 @@
         return;
   }
-  
+
+#if 1
   int tries = min(20, nrays);
 
@@ -1783 +1776 @@
         r1 = Random(nrays);
         r2 = Random(nrays);
-        if (leaf->rays[r1].mRay->mTerminationObject != leaf->rays[r2].mRay->mTerminationObject)
+        if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
           break;
   }
@@ -1789 +1782 @@
   if (r1 == r2)
         r2 = (r1+1)%leaf->rays.size();
-}
-
+
+#else
+  // pick a random ray
+  int base = Random(nrays);
+  RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
+  Intersectable *baseObject = baseRay->GetObject();
+  
+  // and a random 5D derivative which will be used to find the minimal projected distances
+  Vector3 spatialDerivative;
+  Vector3 directionalDerivative;
+
+  while (1) {
+        spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
+                                                                RandomValue(-1.0f, 1.0f),
+                                                                RandomValue(-1.0f, 1.0f));
+        float m = Magnitude(spatialDerivative);
+        if (m != 0) {
+          spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
+          break;
+        }
+  }
+
+  while (1) {
+        directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
+                                                                        RandomValue(-1.0f, 1.0f),
+                                                                        0.0f);
+        float m = Magnitude(directionalDerivative);
+        if (m != 0) {
+          directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
+          break;
+        }
+  }
+
+  // now find the furthest sample from the same object and the closest from a different object
+  int i = 0;
+  float minDist = MAX_FLOAT;
+  float maxDist = -MAX_FLOAT;
+  r1 = base;
+  r2 = base;
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+          ri != leaf->rays.end();
+          ri++, i++) {
+        float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
+                                                                                                         *ri,
+                                                                                                         spatialDerivative,
+                                                                                                         directionalDerivative
+                                                                                                         );
+
+        if ((*ri).GetObject() == baseObject) {
+          if (dist > maxDist) {
+                maxDist = dist;
+                r1 = i;
+          }
+        } else {
+          if (dist > 0.0f && dist < minDist) {
+                minDist = dist;
+                r2 = i;
+          }
+        }
+  }
+  
+#endif
+}
+
+
+struct RayNeighbor {
+  int rayInfo;
+  float distance;
+  RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
+};
+
+void
+RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
+                                                        vector<RssTreeLeaf::SilhouetteRays> &rays
+                                                        )
+{
+  // for every leaf find its neares neighbor from a different object
+  vector<RayNeighbor> neighbors(leaf->rays.size());
+  int i, j;
+  for (i=0; i < leaf->rays.size(); i++)
+        for (j=0; j < leaf->rays.size(); j++)
+          if (i!=j) {
+                float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
+                if (d < neighbors[i].distance) {
+                  neighbors[i].distance = d;
+                  neighbors[i].rayInfo = j;
+                }
+          }
+
+  // now check which are the pairs of nearest neighbors
+  for (i=0; i < leaf->rays.size(); i++) {
+        int j = neighbors[i].rayInfo;
+        if (neighbors[j].rayInfo == i) {
+          // this is a silhouette edge pair
+          if (i < j) {
+                // generate an silhoutte ray pair
+                rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
+                                                                                                   &leaf->rays[j]));
+          }
+        } else {
+          // this is not a silhouette ray - delete???
+        }
+  }
+  
+}
+
+
+bool
+GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
+                                 const int index,
+                                 int &i1,
+                                 int &i2,
+                                 int &i3)
+{
+  int found = 0;
+  int indices[3];
+
+  int size = rays.size();
+  // use russian roulete selection for the tripple
+  // number of free positions for the bullet
+  int positions = size - 1;
+  int i;
+  for (i=0; i < size; i++) {
+        if (rays[i]->mRay->Mailed())
+          positions--;
+  }
+  
+  if (positions < 3)
+        return false;
+  
+  for (i=0; i < size; i++) {
+        if (i != index && !rays[i]->mRay->Mailed()) {
+          float p = (3 - found)/(float)positions;
+          if (Random(1.0f) < p) {
+                indices[found] = i;
+                found++;
+          }
+        }
+        positions--;
+  }
+}
+
+bool
+RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
+                                                                const RssTreeNode::RayInfo &r1,
+                                                                const RssTreeNode::RayInfo &r2,
+                                                                const RssTreeNode::RayInfo &r3)
+{
+  
+
+  return false;
+}
+
+int
+RssTree::RemoveInteriorRays(
+                                                        RssTreeLeaf *leaf
+                                                        )
+{
+#if 1
+  // first collect all objects refered in this leaf
+  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> > rayMap;
+  
+  RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
+  for (; it != leaf->rays.end(); ++it) {
+        Intersectable *object = (*it).GetObject();
+        
+        rayMap[object].push_back(&(*it));
+//      vector<RayInfo *> *data = rayMap.Find(object);
+//      if (data) {
+//        data->push_back(&(*it));
+//      } else {
+//        //      rayMap[object] = vector<RayInfo *>;
+//        rayMap[object].push_back(&(*it));
+//      }
+  }
+
+  // now go through all objects
+  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> >::iterator mi;
+
+  // infos of mailed rays are scheduled for removal
+  VssRay::NewMail();
+  for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
+        vector<RssTreeLeaf::RayInfo *> &rays = (*mi).second;
+        vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
+        int rayIndex = 0;
+        
+        for (; ri != rays.end(); ++ri, ++rayIndex) {
+          RssTreeNode::RayInfo *ray = *ri;
+          int tries = rays.size();
+          for (int i = 0; i < tries; i++) {
+                int r1, r2, r3;
+                GetRandomTripple(rays,
+                                                 rayIndex,
+                                                 r1,
+                                                 r2,
+                                                 r3);
+                if (IsRayConvexCombination(*ray,
+                                                                   *rays[r1],
+                                                                   *rays[r2],
+                                                                   *rays[r3])) {
+                  ray->mRay->Mail();
+                }
+          }
+        }
+  }
+  
+
+#endif  
+        return 0;
+}
 
 void
@@ -1798 +1999 @@
 {
 
+  if (numberOfRays == 0)
+        return;
+  
   int nrays = leaf->rays.size();
-  for (int i=0; i < numberOfRays; i++) {
-        bool useExtendedConvexCombination = (nrays >= 2) && (i > numberOfRays/2);
-
-                
+
+
+  AxisAlignedBox3 box = GetBBox(leaf);
+  AxisAlignedBox3 dirBox = GetDirBBox(leaf);
+
+  AxisAlignedBox3 sBox(box);
+  AxisAlignedBox3 sDirBox(dirBox);
+  const float smoothRange = 0.5f;
+  sBox.Scale(1.0f + smoothRange);
+  sDirBox.Scale(1.0f + smoothRange);
+  int smoothRays = (int)numberOfRays*0.0f;
+  
+#if AVG_LEN
+  float avgLength = 0.0f;
+
+  if (nrays) {
+        const int numAvg = 5;
+        
+        int step = (nrays-1)/10;
+        if (step<1)
+          step = 1;
+        
+        int summed = 0;
+        float sumLength = 0.0f;
+        for (int i=0; i < nrays; i+=step) {
+          sumLength += leaf->rays[i].mRay->GetSize();
+          summed++;
+        }
+  
+        avgLength = sumLength/summed;
+  }
+#endif
+  
+  Exporter *exporter = NULL;
+  VssRayContainer selectedRays;
+  static int counter = 0;
+  if (counter < 0) {
+        char filename[128];
+        sprintf(filename, "selected-rays%04d.x3d", counter);
+        exporter = Exporter::GetExporter(filename);
+        //      exporter->SetWireframe();
+        //      exporter->ExportKdTree(*mKdTree);
+        exporter->SetWireframe();
+        exporter->ExportScene(preprocessor->mSceneGraph->mRoot);
+        exporter->SetWireframe();
+        exporter->ExportBox(box);
+        exporter->SetFilled();
+        counter++;
+  }
+  
+  int numberOfTries = numberOfRays*4;
+  int generated = 0;
+
+  Intersectable::NewMail();
+  vector<Intersectable *> pvs(0);
+  pvs.reserve(leaf->GetPvsSize());
+  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+          ri != leaf->rays.end();
+          ri++) {
+        Intersectable *object = (*ri).GetObject();
+        if (object && !object->Mailed()) {
+          object->Mail();
+          pvs.push_back(object);
+          if (exporter)
+                exporter->ExportIntersectable(object);
+        }
+  }
+  
+  for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
+        bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) < 0.7f));
+        
         Vector3 origin, direction;
         // generate half of convex combination and half of random rays
         if (useExtendedConvexCombination) {
-          // pickup 3 random rays
+          // pickup 2 random rays
           int r1, r2;
           PickEdgeRays(leaf, r1, r2);
+
           
           Vector3 o1 = leaf->rays[r1].GetOrigin();
           Vector3 o2 = leaf->rays[r2].GetOrigin();
           
-          const float overlap = 0.1f;
+          const float overlap = 0.0f;
           float w1, w2;
           GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
           origin = w1*o1 + w2*o2;
+          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
           direction =
-                w1*leaf->rays[r1].mRay->GetDir() +
-                w2*leaf->rays[r2].mRay->GetDir();
-          // shift the origin a little bit
-          origin += direction*0.1f;
+                w1*leaf->rays[r1].GetDir() +
+                w2*leaf->rays[r2].GetDir();
         } else {
-          origin = GetBBox(leaf).GetRandomPoint();
-          Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
+          Vector3 dirVector;
+
+          Vector3 pVector;
+          Vector3 dVector;
+
+          bool useHalton = true;
+          
+          if (useHalton) {
+                // generate a random 5D vector
+                pVector = Vector3(leaf->halton.GetNumber(1),
+                                                  leaf->halton.GetNumber(2),
+                                                  leaf->halton.GetNumber(3));
+                
+                dVector = Vector3(leaf->halton.GetNumber(4),
+                                                  leaf->halton.GetNumber(5),
+                                                  0.0f);
+                leaf->halton.GenerateNext();
+          } else {
+                pVector = Vector3(RandomValue(0.0f, 1.0f),
+                                                  RandomValue(0.0f, 1.0f),
+                                                  RandomValue(0.0f, 1.0f));
+                
+                dVector = Vector3(RandomValue(0.0f, 1.0f),
+                                                  RandomValue(0.0f, 1.0f),
+                                                  0.0f);
+          }
+          
+          if (i < smoothRays) {
+                origin = sBox.GetPoint(pVector);
+                dirVector = sDirBox.GetPoint(dVector);
+          } else {
+                origin = box.GetPoint(pVector);
+                dirVector = dirBox.GetPoint(dVector);
+          }
           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));
-  }
+
+        // shift the origin a little bit
+        direction.Normalize();
+        
+        //      float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
+        float dist = 0.0f;
+        float minT, maxT;
+        // compute interection of the ray with the box
+
+        if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
+          dist = maxT;
+  
+        
+        origin += direction*dist;
+        
+        bool intersects = false;
+        
+
+        if (i > numberOfRays/2) {
+          if (exporter) {
+                VssRay *ray = new VssRay(origin, origin + 100*direction, NULL, NULL);
+                selectedRays.push_back(ray);
+          }
+
+          // check whether the ray does not intersect already visible objects
+          Ray traversalRay;
+          traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
+          for(vector<Intersectable *>::const_iterator oi = pvs.begin();
+                  oi != pvs.end();
+                  oi++) {
+                Intersectable *object = *oi;
+                
+                if ( object->CastRay(traversalRay) ) {
+                  intersects = true;
+                  break;
+                }
+          }
+        }
+        
+        
+        if (!intersects) {
+          //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
+          rays.push_back(SimpleRay(origin, direction));
+          generated++;
+        }
+  }
+
+  //  cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
+  
+  if (exporter) {
+        exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
+        delete exporter;
+        CLEAR_CONTAINER(selectedRays);
+  }
+}
+
+int
+RssTree::PruneRays(RssTreeLeaf *leaf,
+                                   const float contributionThreshold)
+{
+  int i;
+  int j;
+  
+  for (j=0, i=0; i < leaf->rays.size(); i++) {
+        
+        if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
+          // copy a valid sample
+          if (i!=j) 
+                leaf->rays[j] = leaf->rays[i];
+          j++;
+        } else {
+          // delete the ray
+          leaf->rays[i].mRay->Unref();
+          if (leaf->rays[i].mRay->RefCount() != 0) {
+                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
+                exit(1);
+          }
+          delete leaf->rays[i].mRay;
+        }
+  }
+
+
+  leaf->rays.resize(j);
+  int removed = (i-j);
+  stat.rayRefs -= removed;
+  return removed;
+}
+
+int
+RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
+                                                 const float ratio)
+{
+  int i;
+  int j;
+  
+  for (j=0, i=0; i < leaf->rays.size(); i++) {
+        
+        if (Random(1.0f) < ratio) {
+          // copy a valid sample
+          if (i!=j) 
+                leaf->rays[j] = leaf->rays[i];
+          j++;
+        } else {
+          // delete the ray
+          leaf->rays[i].mRay->Unref();
+          if (leaf->rays[i].mRay->RefCount() != 0) {
+                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
+                exit(1);
+          }
+          delete leaf->rays[i].mRay;
+        }
+  }
+
+
+  leaf->rays.resize(j);
+  int removed = (i-j);
+  stat.rayRefs -= removed;
+  return removed;
+}
+
+int
+RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
+                                                           const float ratio)
+{
+  int i;
+  
+  if (leaf->rays.size() == 0)
+        return 0;
+  
+  sort(leaf->rays.begin(),
+           leaf->rays.end(),
+           RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
+
+  int desired = ratio*leaf->rays.size();
+  int removed = leaf->rays.size() - desired;
+  
+  for (i=desired; i < leaf->rays.size(); i++) {
+        // delete the ray
+        leaf->rays[i].mRay->Unref();
+        if (leaf->rays[i].mRay->RefCount() != 0) {
+          cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
+          exit(1);
+        }
+        delete leaf->rays[i].mRay;
+  }
+  
+  leaf->rays.resize(desired);
+  stat.rayRefs -= removed;
+  return removed;
+}
+
+
+int
+RssTree::PruneRays(
+                                   const int desired
+                                   )
+{
+  bool globalPrunning = false;
+
+  stack<RssTreeNode *> tstack;
+  int prunned = 0;
+
+  Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl;
+
+  if (globalPrunning) {
+        VssRayContainer allRays;
+        allRays.reserve(stat.rayRefs);
+        CollectRays(allRays);
+        sort(allRays.begin(),
+                 allRays.end(),
+                 GreaterWeightedPvsContribution);
+        
+        if ( desired >= allRays.size() )
+          return 0;
+        
+        float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
+        
+        tstack.push(root);
+        
+        while (!tstack.empty()) {
+          RssTreeNode *node = tstack.top();
+          tstack.pop();
+          
+          if (node->IsLeaf()) {
+                RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+                prunned += PruneRays(leaf, contributionThreshold);
+          } else {
+                RssTreeInterior *in = (RssTreeInterior *)node;
+                // both nodes for directional splits
+                tstack.push(in->front);
+                tstack.push(in->back);
+          }
+        }
+  } else {
+        // prune random rays from each leaf so that the ratio's remain the same
+        tstack.push(root);
+        float ratio = desired/(float)stat.rayRefs;
+        
+        while (!tstack.empty()) {
+          RssTreeNode *node = tstack.top();
+          tstack.pop();
+          
+          if (node->IsLeaf()) {
+                RssTreeLeaf *leaf = (RssTreeLeaf *)node;
+                //              prunned += PruneRaysRandom(leaf, ratio);
+                prunned += PruneRaysContribution(leaf, ratio);
+          } else {
+                RssTreeInterior *in = (RssTreeInterior *)node;
+                // both nodes for directional splits
+                tstack.push(in->front);
+                tstack.push(in->back);
+          }
+        }
+
+        
+
+
+  }
+  
+  
+  
+  Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl;
+  
+  return prunned;
 }
 
@@ -1837 +2361 @@
                                           SimpleRayContainer &rays)
 {
+
 
   vector<RssTreeLeaf *> leaves;
@@ -1859 +2384 @@
                                 
   float avgContrib = sumContrib/numberOfLeaves;
-  float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
+  
+  // always generate at leat n ray per leaf
+  int fixedPerLeaf = 1;
+  int fixed = 1*leaves.size();
+  int iGenerated = numberOfRays;
+  float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
+
   k = 0;
   for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
@@ -1866 +2397 @@
           RssTreeLeaf *leaf = leaves[i];
           float c = leaf->GetImportance();
-          int num = (c*ratioPerLeaf + 0.5);
+
+          int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
           GenerateLeafRays(leaf, num, rays);
         }
@@ -1889 +2421 @@
           RssTreeLeaf *leaf = (RssTreeLeaf *)node;
           // update pvs size
-          leaf->UpdatePvsSize();
+          UpdatePvsSize(leaf);
           sumPvs += leaf->GetPvsSize();
           leaves++;
@@ -1904 +2436 @@
 }
 
-        
+float weightAbsContributions = 0.0f;
+// if small very high importance of the last sample
+// if 1.0f then weighs = 1 1/2 1/3 1/4
+float passSampleWeightDecay = 0.5f;
+
+void
+RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
+                                                        float &weight,
+                                                        float &contribution)
+{
+  VssRay *ray = info.mRay;
+  weight = 1.0f/(mCurrentPass - ray->mPass + passSampleWeightDecay);
+  contribution =
+        weightAbsContributions*ray->mPvsContribution +
+        (1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
+  // store the computed value
+  info.mRay->mWeightedPvsContribution = weight*contribution;
+}
+
+float
+RssTree::GetSampleWeight(const int pass)
+{
+  int passDiff = mCurrentPass - pass;
+  float weight;
+  if (1) 
+        weight = 1.0f/(passDiff + passSampleWeightDecay);
+  else 
+        switch (passDiff) {
+        case 0:
+          weight = 1.0f;
+          break;
+        default:
+          weight = 0.01f;
+          break;
+          //    case 1:
+//        weight = 0.5f;
+//        break;
+//      case 2:
+//        weight = 0.25f;
+//        break;
+//      case 3:
+//        weight = 0.12f;
+//        break;
+//      case 4:
+//        weight = 0.06f;
+//        break;
+//      default:
+//        weight = 0.03f;
+//        break;
+        }
+  return weight;
+}
+
+void
+RssTree::ComputeImportance(RssTreeLeaf *leaf) 
+{
+  if (0) 
+        leaf->mImportance = leaf->GetAvgRayContribution();
+  else {
+        RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
+          it_end = leaf->rays.end();
+        
+        float sumContributions = 0.0f;
+        float sumRelContributions = 0.0f;
+        float sumWeights = 0.0f;
+        for (; it != it_end; ++it) {
+          VssRay *ray = (*it).mRay;
+          float weight = GetSampleWeight(ray->mPass);
+          
+          sumContributions += weight*ray->mPvsContribution;
+          sumRelContributions += weight*ray->mRelativePvsContribution;
+          
+          //      sumWeights += weight;
+          sumWeights += 1.0f;
+        }
+        // $$ 
+        // sumWeights = leaf->mTotalRays;
+          
+        if (sumWeights != 0.0f) 
+          leaf->mImportance =
+                (weightAbsContributions*sumContributions +
+                 (1.0f - weightAbsContributions)*sumRelContributions)/sumWeights;
+        else
+          leaf->mImportance = 0.0f;
+        
+  }
+  
+  // return GetAvgRayContribution()*mImportance;
+  //return GetAvgRayContribution();
+}
+
+
 float
 RssTreeLeaf::GetImportance()  const
 {
-
-  if (1) {
-        return GetAvgRayContribution();
-        //      return GetPvsSize();
-  } else {
-        // return GetAvgRayContribution()*mEntropyImportance;
-        //return GetAvgRayContribution();
-        return mEntropyImportance;
-  }
+  return mImportance;
 }
 
 
 float
-RssTreeLeaf::ComputePvsEntropy()
+RssTreeLeaf::ComputePvsEntropy() const
 {
   int samples = 0;
   Intersectable::NewMail();
   // set all object as belonging to the fron pvs
-  for(RssTreeNode::RayInfoContainer::iterator ri = rays.begin();
+  for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
           ri != rays.end();
           ri++)
         if ((*ri).mRay->IsActive()) {
-          Intersectable *object = (*ri).mRay->mTerminationObject;
+          Intersectable *object = (*ri).GetObject();
           if (object) {
                 if (!object->Mailed()) {
@@ -1949 +2564 @@
                 ri++) 
           if ((*ri).mRay->IsActive()) {
-                Intersectable *object = (*ri).mRay->mTerminationObject;
+                Intersectable *object = (*ri).GetObject();
                 if (object) {
                   if (!object->Mailed()) {
@@ -1967 +2582 @@
 
 float
-RssTreeLeaf::ComputeRayLengthEntropy()
+RssTreeLeaf::ComputeRayLengthEntropy() const
 {
   // get sum of all ray lengths
@@ -2017 +2632 @@
 
 
-void
-RssTreeLeaf::ComputeEntropyImportance()
-{
-  mPvsEntropy = ComputePvsEntropy();
-  mRayLengthEntropy = ComputeRayLengthEntropy();
+float
+RssTreeLeaf::ComputeEntropyImportance() const
+{
   
   //  mEntropy = 1.0f - ComputeRayLengthEntropy();
-  mEntropyImportance = 1.0f - ComputePvsEntropy();
-  
-  //  cout<<"ei="<<mEntropyImportance<<" ";
+  return 1.0f - ComputePvsEntropy();
+}
+
+float
+RssTreeLeaf::ComputeRayContributionImportance() const
+{
+  //  mPvsEntropy = ComputePvsEntropy();
+  //  mRayLengthEntropy = ComputeRayLengthEntropy();
+  
+  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
+  return 1.0f - ComputePvsEntropy();
 }
 
@@ -2059 +2680 @@
   VssRayContainer allRays;
   CollectRays(allRays);
+
   
   int desired = min(number, (int)allRays.size());
@@ -2105 +2727 @@
   return (int)rays.size();
 }
+

trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.h

@@ -21 +21 @@
 #include "VssRay.h"
 #include "AxisAlignedBox3.h"
+#include "Halton.h"
 
 
@@ -69 +70 @@
   // number of dynamically removed ray refs
   int removedRayRefs;
-  
+
+
+
+  // dynamically updated values
+  float avgPvsSize;
+  float avgRays;
+  float avgRayContribution;
+  float avgPvsEntropy;
+  float avgRayLengthEntropy;
+  float avgImportance;
+
+  float avgWeightedRayContribution;
+
   // Constructor
   RssStatistics() {
@@ -144 +157 @@
           return a.mRay < b.mRay;
         }
+
+#define USE_ORIGIN 0
         
         Vector3 GetOrigin() const {
+#if USE_ORIGIN
           return mRay->GetOrigin();
-        }
-
-        float ExtrapOrigin(const int axis) const {
+#else
+          return mRay->GetTermination();
+#endif
+        }
+
+        Vector3 GetTermination() const {
+#if USE_ORIGIN
+          return mRay->GetTermination();
+#else
+          return mRay->GetOrigin();
+#endif
+        }
+
+        float GetOrigin(const int axis) const {
+#if USE_ORIGIN
           return mRay->GetOrigin(axis);
-        }
-        
-        float ExtrapTermination(const int axis) const {
+#else
           return mRay->GetTermination(axis);
-        }
-
-        Vector3 Extrap(const float t) const {
-          return mRay->Extrap(t);
-        }
+#endif
+        }
+
+        Vector3 GetDir() const {
+#if USE_ORIGIN
+          return mRay->GetDir();
+#else
+          return -mRay->GetDir();
+#endif
+        }
+
+        float GetDirParametrization(const int axis) const {
+#if USE_ORIGIN
+          return mRay->GetDirParametrization(axis);
+#else
+          return mRay->GetOpositeDirParametrization(axis);
+#endif
+        }
+
+        
+//      Vector3 GetTermination() const {
+//        return mRay->GetOrigin();
+//      }
+        
+//      float GetTermination(const int axis) const {
+//        return mRay->GetTermination(axis);
+//      }
+
+        Intersectable *GetObject() const {
+#if USE_ORIGIN
+          return mRay->mTerminationObject;
+#else
+          return mRay->mOriginObject;
+#endif
+        }
+        
+        Intersectable *GetSourceObject() const {
+#if USE_ORIGIN
+          return mRay->mOriginObject;
+#else
+          return mRay->mTerminationObject;
+#endif
+        }
+
+        //      Vector3 Extrap(const float t) const {
+        //        return mRay->Extrap(t);
+        //      }
         
         int
@@ -167 +235 @@
           
           // only compare the side of the origin
-          float rpos = mRay->GetOrigin(axis);
+          float rpos = (axis < 3) ? GetOrigin(axis) : GetDirParametrization(axis - 3);
           return (rpos <= position) ? -1 : 1;
         }
-        
+
+        static bool GreaterWeightedPvsContribution(const RayInfo &a,
+                                                                                           const RayInfo &b) {
+          return a.mRay->mWeightedPvsContribution > b.mRay->mWeightedPvsContribution;
+        }
+
+        static float Distance(const RayInfo &a,
+                                                 const RayInfo &b)
+        {
+          if (a.mRay->mTerminationObject == b.mRay->mTerminationObject)
+                return MAX_FLOAT;
+          
+          return Min(a.mRay->SqrDistance(b.GetTermination()),
+                                 b.mRay->SqrDistance(a.GetTermination()));
+        }
+
+        
+        static float SqrDistance5D(const RayInfo &a,
+                                                           const RayInfo &b,
+                                                           const Vector3 &spatialDerivative,
+                                                           const Vector3 &directionalDerivative
+                                                           )
+        {
+          Vector3 spatialDiff = b.GetOrigin() - a.GetOrigin();
+          Vector3 directionalDiff = b.GetDir() - a.GetDir();
+          
+          return DotProd(spatialDiff, spatialDerivative) +
+                DotProd(directionalDiff, directionalDerivative);
+        }
+
+        
+  };
+
+  struct SilhouetteRays {
+        RayInfo *a;
+        RayInfo *b;
+        SilhouetteRays() {}
+        SilhouetteRays(RayInfo *_a, RayInfo *_b):a(_a), b(_b) {
+        }
   };
   
@@ -307 +413 @@
   
   RayInfoContainer rays;
-  int mPassingRays;
+  int mTotalRays;
         
   bool mValidPvs;
-  float mEntropyImportance;
-  float mPvsEntropy;
-  float mRayLengthEntropy;
-  
-        
+  float mImportance;
+  
+  HaltonSequence halton;
+  
   RssTreeLeaf(RssTreeInterior *p,
                           const int nRays
-                          ):RssTreeNode(p), rays(), mPvsSize(0), mPassingRays(0), mValidPvs(false) {
+                          ):RssTreeNode(p), rays(), mPvsSize(0), mTotalRays(0), mValidPvs(false) {
     rays.reserve(nRays);
   }
@@ -328 +433 @@
     s<<
           "L: rays="<<(int)rays.size()<<
-          " pvs="<<mPvsSize<<" pEntropy="<<mPvsEntropy<<" lEntropy="<<mRayLengthEntropy<<endl;
+          " pvs="<<mPvsSize<<" importance="<<mImportance<<endl;
   };
   
@@ -334 +439 @@
         mValidPvs = false;
     rays.push_back(data);
-    data.mRay->Ref();
-        if (data.GetRayClass() == RayInfo::PASSING_RAY)
-          mPassingRays++;
+        mTotalRays++;
+        data.mRay->Ref();
   }
         
@@ -348 +452 @@
   }
 
-  void
-  UpdatePvsSize();
 
   float
-  ComputePvsEntropy();
+  ComputePvsEntropy() const;
   
   float
-  ComputeRayLengthEntropy();
+  ComputeRayLengthEntropy() const;
   
   float
-  ComputeRayTerminationEntropy();
-  
-  void
-  ComputeEntropyImportance();
+  ComputeRayTerminationEntropy() const;
+  
+
+  float
+  ComputeEntropyImportance() const;
+  
+  float
+  ComputeRayContributionImportance() const;
 
   void Mail() { mailbox = mailID; }
@@ -401 +507 @@
         return a->GetAvgRayContribution() > b->GetAvgRayContribution();
   }
+
   
 };
@@ -560 +667 @@
   // depth at which directional splits are performed if mInterleaveDirSplits is false
   int mDirSplitDepth;
+
+  // maximal number of rays maintained in the rss tree
+  int mMaxRays;
   
   // maximal size of the box on which the refdir splitting can be performed
@@ -582 +692 @@
   int minCollapseDepth;
 
+  bool mImportanceBasedCost;
+
+  // sampling pass
+  int mCurrentPass;
   
   // reusable array of split candidates
@@ -627 +741 @@
   Subdivide(const TraversalData &tdata);
         
-  int
-  SelectPlane(RssTreeLeaf *leaf,
-                          const AxisAlignedBox3 &box,
-                          float &position,
-                          int &raysBack,
-                          int &raysFront,
-                          int &pvsBack,
-                          int &pvsFront
-                          );
 
   void
@@ -658 +763 @@
   }
   
-
-  float
+  struct SplitInfo {
+        
+        int axis;
+        float position;
+        float costRatio;
+        
+        int rays;
+        int raysBack;
+        int raysFront;
+
+        int pvs;
+        int pvsBack;
+        int pvsFront;
+
+        float contribution;
+        float contributionBack;
+        float contributionFront;
+
+        int viewCells;
+        int viewCellsBack;
+        int viewCellsFront;
+
+        float volume;
+        float volumeBack;
+        float volumeFront;
+        
+  };
+  
+  void
   BestCostRatio(
                                 RssTreeLeaf *node,
-                                int &axis,
-                                float &position,
-                                int &raysBack,
-                                int &raysFront,
-                                int &pvsBack,
-                                int &pvsFront
+                                SplitInfo &info
+//                              int &axis,
+//                              float &position,
+//                              int &raysBack,
+//                              int &raysFront,
+//                              int &pvsBack,
+//                              int &pvsFront
                                 );
         
-  float
+  
+  void
   EvalCostRatio(
                                 RssTreeLeaf *node,
-                                const int axis,
-                                const float position,
-                                int &raysBack,
-                                int &raysFront,
-                                int &pvsBack,
-                                int &pvsFront
+                                SplitInfo &info
+//                              const int axis,
+//                              const float position,
+//                              int &raysBack,
+//                              int &raysFront,
+//                              int &pvsBack,
+//                              int &pvsFront
                                 );
 
-  float
+  void
   EvalCostRatioHeuristic(
                                                  RssTreeLeaf *node,
-                                                 const int axis,
-                                                 float &position,
-                                                 int &raysBack,
-                                                 int &raysFront,
-                                                 int &pvsBack,
-                                                 int &pvsFront
+                                                 SplitInfo &info
+//                                               const int axis,
+//                                               float &position,
+//                                               int &raysBack,
+//                                               int &raysFront,
+//                                               int &pvsBack,
+//                                               int &pvsFront
                                                  );
 
-  float
+  int
+  SelectPlane(RssTreeLeaf *leaf,
+                          const AxisAlignedBox3 &box,
+                          SplitInfo &info
+                          );
+  
+  void
   GetCostRatio(
                            RssTreeLeaf *leaf,
-                           const int axis,
-                           const float position,
-                           const int raysBack,
-                           const int raysFront,
-                           const int pvsBack,
-                           const int pvsFront
+                           SplitInfo &info
+                           //                      const int axis,
+                           //                      const float position,
+                           //                      const int raysBack,
+                           //                      const int raysFront,
+                           //                      const int pvsBack,
+                           //                      const int pvsFront
                            );
-
+  
   AxisAlignedBox3 GetBBox(const RssTreeNode *node) const {
     if (node->parent == NULL)
@@ -803 +946 @@
         
   void
-  GetTreeStatistics(
-                                        float &avgPvsSize,
-                                        float &avgRays,
-                                        float &avgRayContribution,
-                                        float &avgPvsEntropy,
-                                        float &avgRayLengthEntropy,
-                                        float &avgImportance
-                                        );
+  UpdateTreeStatistics();
 
 
@@ -865 +1001 @@
                           );
 
+  void
+  UpdatePvsSize(RssTreeLeaf *leaf);
+
+  void
+  ComputeImportance(RssTreeLeaf *leaf);
+  
+  void
+  SetPass(const int p) {
+        mCurrentPass = p;
+  }
+  
+  void GetRayContribution(const RssTreeNode::RayInfo &info,
+                                                  float &weight,
+                                                  float &contribution);
+
+  float
+  GetSampleWeight(const int pass);
+
+  int
+  RemoveInteriorRays(
+                                         RssTreeLeaf *leaf
+                                         );
+
+  bool
+  IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
+                                                 const RssTreeNode::RayInfo &r1,
+                                                 const RssTreeNode::RayInfo &r2,
+                                                 const RssTreeNode::RayInfo &r3);
+
+  int
+  PruneRays(RssTreeLeaf *leaf,
+                        const float contributionThreshold);
+  int
+  PruneRaysRandom(RssTreeLeaf *leaf,
+                                  const float ratio);
+
+  int
+  PruneRaysContribution(RssTreeLeaf *leaf,
+                                                const float ratio);
+
+  int
+  PruneRays(
+                        const int desired
+                        );
+  
+
+  void
+  FindSilhouetteRays(RssTreeLeaf *leaf,
+                                         vector<RssTreeLeaf::SilhouetteRays> &rays
+                                         );
+  
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp

@@ -242 +242 @@
 SamplingPreprocessor::ComputeVisibility()
 {
-  // pickup an object
-  ObjectContainer objects;
   
   /// rays per pass
   RayContainer passRays;
 
-  mSceneGraph->CollectObjects(&objects);
 
   mViewCellsManager->SetViewSpaceBox(mKdTree->GetBox());
@@ -270 +267 @@
                 int reverseSamples = 0;
                 
-                for (i = 0; i < objects.size(); i++) {
+                for (i = 0; i < mObjects.size(); i++) {
                                                 
                         KdNode *nodeToSample = NULL;
-                        Intersectable *object = objects[i];
+                        Intersectable *object = mObjects[i];
                 
                         int pvsSize = 0;
@@ -339 +336 @@
                 // cast rays for view cells construction
                 ProcessViewCells(passRays,
-                                                 objects,
+                                                 mObjects,
                                                  passSampleContributions,
                                                  passContributingSamples);
@@ -361 +358 @@
                                  << "k   #SampleContributions=" << passSampleContributions << " (" 
                                  << 100*passContributingSamples/(float)passSamples<<"%)" << " avgPVS="
-                                 << (float)pvsSize /(float)objects.size() << endl 
+                                 << (float)pvsSize /(float)mObjects.size() << endl 
                                  << "avg ray contrib=" << avgRayContrib << endl
                                  << "reverse samples [%]" << reverseSamples/(float)passSamples*100.0f << endl;
@@ -371 +368 @@
                         "#SampleContributions\n" << passSampleContributions << endl << 
                         "#PContributingSamples\n"<<100*passContributingSamples/(float)passSamples<<endl <<
-                        "#AvgPVS\n"<< pvsSize/(float)objects.size() << endl <<
+                        "#AvgPVS\n"<< pvsSize/(float)mObjects.size() << endl <<
                         "#AvgRayContrib\n" << avgRayContrib << endl;
         }
@@ -387 +384 @@
         // construct view cells if not already constructed
         if (!mViewCellsManager->ViewCellsConstructed())
-                mViewCellsManager->Construct(objects, mVssSampleRays);
+          mViewCellsManager->Construct(mObjects, mVssSampleRays);
 
         // $$JB temporary removed

trunk/VUT/GtpVisibilityPreprocessor/src/SceneGraph.cpp

@@ -5 +5 @@
 #include "X3dExporter.h"
 #include "Intersectable.h"
-
 
 
@@ -104 +103 @@
         
 }
+
+
+void
+SceneGraphNode::UpdateBox()
+{
+  AxisAlignedBox3 box;
+  
+  box.Initialize();
+  
+  ObjectContainer::const_iterator mi = mGeometry.begin();
+  for (; mi != mGeometry.end(); mi++)
+        box.Include((*mi)->GetBox());
+  
+  SceneGraphNodeContainer::iterator ni = mChildren.begin();
+  for (; ni != mChildren.end(); ni++) {
+        (*ni)->UpdateBox();
+        box.Include((*ni)->mBox);
+  }
+  
+  mBox = box;
+}

trunk/VUT/GtpVisibilityPreprocessor/src/SceneGraph.h

@@ -16 +16 @@
   SceneGraphNodeContainer mChildren;
   AxisAlignedBox3 mBox;
+  void UpdateBox();
 };
 
@@ -27 +28 @@
   
   int
-        CollectObjects(ObjectContainer *instances);
-        
+  CollectObjects(ObjectContainer *instances);
+  
   int
-        AssignObjectIds();
+  AssignObjectIds();
+  
+  void
+  GetStatistics(int &intersectables, int &faces) const;
 
-        void
-        GetStatistics(int &intersectables, int &faces) const;
-
+  AxisAlignedBox3 GetBox() const { return mRoot->mBox; }
+  
 protected:
 };
@@ -40 +43 @@
 
 
-
 #endif

trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp

@@ -191 +191 @@
 
 Vector3
-UniformRandomVector(const Vector3 &normal)
-{
-        //  float r1 = RandomValue(0.0f, 1.0f);
-        //  float r2 = RandomValue(0.0f, 1.0f);
-        float r1, r2;
+UniformRandomVector()
+{
+  //  float r1 = RandomValue(0.0f, 1.0f);
+  //  float r2 = RandomValue(0.0f, 1.0f);
+  float r1, r2;
         
-        halton2.GetNext(r1, r2);
+  halton2.GetNext(r1, r2);
         
         
-        float cosTheta = 1.0f - r1;
+  float cosTheta = 1.0f - 2*r1;
   float sinTheta = sqrt(1 - sqr(cosTheta));
   float fi = 2.0f*M_PI*r2;
   
   Vector3 dir(sinTheta*sin(fi),
-                                                        cosTheta,
-                                                        sinTheta*cos(fi));
-        
+                          cosTheta,
+                          sinTheta*cos(fi));
+  
+  return dir;
+}
+
+Vector3
+UniformRandomVector(const Vector3 &normal)
+{
+  //  float r1 = RandomValue(0.0f, 1.0f);
+  //  float r2 = RandomValue(0.0f, 1.0f);
+  float r1, r2;
+  
+  halton2.GetNext(r1, r2);
+  
+  
+  float cosTheta = 1.0f - r1;
+  float sinTheta = sqrt(1 - sqr(cosTheta));
+  float fi = 2.0f*M_PI*r2;
+  
+  Vector3 dir(sinTheta*sin(fi),
+                          cosTheta,
+                          sinTheta*cos(fi));
+  
   
   //  return Normalize(dir);
-        
+  
   Matrix4x4 m = RotationVectorsMatrix(
-                                                                                                                                                        normal,
-                                                                                                                                                        Vector3(0,1,0));
+                                                                          normal,
+                                                                          Vector3(0,1,0));
   Matrix4x4 mi = Invert(m);
   m = m*RotationVectorsMatrix(
-                                                                                                                        Vector3(0,1,0),
-                                                                                                                        Normalize(dir))*mi;
+                                                          Vector3(0,1,0),
+                                                          Normalize(dir))*mi;
   
   return TransformNormal(m, normal);
   
-        //    return TransformNormal(
-        //                       RotationVectorsMatrix(
-        //                                             Vector3(0,1,0),
-        //                                             Normalize(dir)
-        //                                             ),
-        //                       normal
-        //                       );
-}
+  //    return TransformNormal(
+  //                     RotationVectorsMatrix(
+  //                                           Vector3(0,1,0),
+  //                                           Normalize(dir)
+  //                                           ),
+  //                     normal
+  //                     );
+}
+

trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h

@@ -228 +228 @@
 
   friend Vector3 UniformRandomVector(const Vector3 &normal);
+  friend Vector3 UniformRandomVector();
+  
+
 };
 
@@ -445 +448 @@
 }
 
+
 // Overload << operator for C++-style output
 inline ostream&

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp

@@ -2672 +2672 @@
 }
 
+
+ViewCell *
+BspTree::GetViewCell(const Vector3 &point)
+{
+  stack<BspNode *> nodeStack;
+  nodeStack.push(mRoot);
+  
+  ViewCell *viewcell = NULL;
+  
+  while (!nodeStack.empty())  {
+        BspNode *node = nodeStack.top();
+        nodeStack.pop();
+        
+        if (node->IsLeaf()) {
+          viewcell = dynamic_cast<BspLeaf *>(node)->mViewCell;
+          break;
+        } else {
+          
+          BspInterior *interior = dynamic_cast<BspInterior *>(node);
+                
+          // random decision
+          if (interior->GetPlane().Side(point) < 0)
+                nodeStack.push(interior->GetBack());
+          else
+                nodeStack.push(interior->GetFront());
+        }
+  }
+  
+  return viewcell;
+}
+
 void BspNodeGeometry::ComputeBoundingBox(AxisAlignedBox3 &box)
 {

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h

@@ -468 +468 @@
                                   );
 
+  ViewCell *
+  GetViewCell(const Vector3 &point);
+  
         /// bsp tree construction types
         enum {FROM_INPUT_VIEW_CELLS, FROM_SCENE_GEOMETRY, FROM_SAMPLES};

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.cpp

@@ -330 +330 @@
 void ViewCellsManager::ComputeSampleContributions(VssRay &ray)
 {
-        ViewCellContainer viewcells;
-                
-        ray.mPvsContribution = 0;
-        ray.mRelativePvsContribution = 0.0f;
-
-        // matt TODO: remove this!!
-        Ray hray(ray);
-        float tmin = 0, tmax = 1.0;
-        
-        //hray.Init(ray.GetOrigin(), ray.GetDir(), Ray::LINE_SEGMENT);
-        if (!GetSceneBbox().GetRaySegment(hray, tmin, tmax) || (tmin > tmax))
-                return;
-
-        Vector3 origin = hray.Extrap(tmin);
-        Vector3 termination = hray.Extrap(tmax);
-        
-        CastLineSegment(origin,
-                                        termination,
-                                        viewcells);
+
+
+
+  ViewCellContainer viewcells;
+                
+  ray.mPvsContribution = 0;
+  ray.mRelativePvsContribution = 0.0f;
+
+  // matt TODO: remove this!!
+  Ray hray(ray);
+  float tmin = 0, tmax = 1.0;
   
-        ViewCellContainer::const_iterator it = viewcells.begin();
-        
-        for (; it != viewcells.end(); ++it) 
-        {       
-                ViewCell *viewcell = *it;
-        
-                // if ray not outside of view space
-                float contribution;
-                bool added = 
-                        viewcell->GetPvs().AddSample(ray.mTerminationObject,
-                                                                                 contribution);
-                
-                if (added)
-                        ++ ray.mPvsContribution;
-                
-                ray.mRelativePvsContribution += contribution;
-        }
+  //hray.Init(ray.GetOrigin(), ray.GetDir(), Ray::LINE_SEGMENT);
+  if (!GetSceneBbox().GetRaySegment(hray, tmin, tmax) || (tmin > tmax))
+        return;
+  
+  Vector3 origin = hray.Extrap(tmin);
+  Vector3 termination = hray.Extrap(tmax);
+  
+  CastLineSegment(origin,
+                                  termination,
+                                  viewcells);
+  
+  // copy viewcells memory efficiently
+  const bool storeViewcells = false;
+  if (storeViewcells) {
+        ray.mViewCells.reserve(viewcells.size());
+        ray.mViewCells = viewcells;
+  }
+  
+  ViewCellContainer::const_iterator it = viewcells.begin();
+
+  bool addInPlace = false;
+
+  if (addInPlace) {
+        for (; it != viewcells.end(); ++it) {   
+          ViewCell *viewcell = *it;
+          
+          // if ray not outside of view space
+          float contribution;
+          bool added = 
+                viewcell->GetPvs().AddSample(ray.mTerminationObject,
+                                                                         contribution
+                                                                         );
+          if (added)
+                ray.mPvsContribution++;
+          
+          ray.mRelativePvsContribution += contribution;
+        }
+  } else {
+    for (; it != viewcells.end(); ++it) {       
+          ViewCell *viewcell = *it;
+          // if ray not outside of view space
+          float contribution;
+          if (viewcell->GetPvs().GetSampleContribution(ray.mTerminationObject,
+                                                                                                   contribution
+                                                                                                   ))
+                ray.mPvsContribution++;
+          
+          ray.mRelativePvsContribution += contribution;
+        }
+
+        for (it = viewcells.begin(); it != viewcells.end(); ++it) {     
+          ViewCell *viewcell = *it;
+          // if ray not outside of view space
+          viewcell->GetPvs().AddSample(ray.mTerminationObject);
+        }
+  }
 }
 
@@ -520 +551 @@
 #endif
 }
-
+  
 
 float BspViewCellsManager::GetRendercost(ViewCell *viewCell, float objRendercost) const
@@ -1705 +1736 @@
                                                                           const VssRayContainer &rays)
 {
-        // if view cells were already constructed
-        if (ViewCellsConstructed())
-                return 0;
-
-        int sampleContributions = 0;
+  // if view cells were already constructed
+  if (ViewCellsConstructed())
+        return 0;
+  
+  Debug << "Constructing bsp view cells" << endl;
+  
+  int sampleContributions = 0;
+  
+  VssRayContainer sampleRays;
+  
+  int limit = min (mConstructionSamples, (int)rays.size());
+
         
         VssRayContainer constructionRays;
@@ -2268 +2306 @@
                 }
         }
-
         return maxDist;
 }
+
+
+
+
+ViewCell *
+VspBspViewCellsManager::GetViewCell(const Vector3 &point)
+{
+  if (!mVspBspTree)
+        return NULL;
+  return mVspBspTree->GetViewCell(point);
+}

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.h

@@ -171 +171 @@
         virtual void GetPvsStatistics(PvsStatistics &stat);
 
-        virtual void PrintPvsStatistics(ostream &s);
+  /** Get a viewcell containing the specified point */
+  virtual ViewCell *GetViewCell(const Vector3 &point) = 0;
+  
+  virtual void
+  PrintPvsStatistics(ostream &s);
 
         /** Returns probability that view point lies in one view cell.
@@ -339 +343 @@
         AxisAlignedBox3 GetSceneBbox() const;
 
+  /** Get a viewcell containing the specified point */
+  ViewCell *GetViewCell(const Vector3 &point) { return NULL; }
+
 protected:
 
@@ -404 +411 @@
         /** Prints out statistics of this approach.
         */
-        //void PrintStatistics(ostream &s) const;
-
-        float GetProbability(ViewCell *viewCell);
-        float GetRendercost(ViewCell *viewCell, float objRendercost) const;
-        
-        AxisAlignedBox3 GetSceneBbox() const;
+  //  virtual void PrintStatistics(ostream &s) const;
+  ViewCell *GetViewCell(const Vector3 &point) { return NULL; }
+
+  float GetProbability(ViewCell *viewCell);
+  float GetRendercost(ViewCell *viewCell, float objRendercost) const;
+  
+  AxisAlignedBox3 GetSceneBbox() const;
 
 protected:
@@ -460 +468 @@
                                                 ViewCellContainer &viewcells);
 
+  ViewCell *GetViewCell(const Vector3 &point) { return NULL; }
+
         float GetProbability(ViewCell *viewCell);
         float GetRendercost(ViewCell *viewCell, float objRendercost) const;
 
         AxisAlignedBox3 GetSceneBbox() const;
-
 protected:
 
@@ -508 +517 @@
         bool ViewCellsConstructed() const;
 
-        //void PrintStatistics(ostream &s) const;
         
         int CastLineSegment(const Vector3 &origin,
@@ -519 +527 @@
         AxisAlignedBox3 GetSceneBbox() const;
 
+  ViewCell *GetViewCell(const Vector3 &point);
+
         bool GetViewPoint(Vector3 &viewPoint) const;
 
         bool ViewPointValid(const Vector3 &viewPoint) const;
-
 protected:
         /** DEPRECATED

trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp

@@ -1525 +1525 @@
 void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
 {
-        stack<BspNode *> nodeStack;
-        nodeStack.push(mRoot);
-
-        ViewCell::NewMail();
-
-        while (!nodeStack.empty())
-        {
-                BspNode *node = nodeStack.top();
-                nodeStack.pop();
-
-                if (node->IsLeaf())
-                {
-                        ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
-                        if (node->TreeValid() && !viewCell->Mailed())
+  stack<BspNode *> nodeStack;
+
+  if (!mRoot)
+        return;
+
+  nodeStack.push(mRoot);
+  
+  ViewCell::NewMail();
+  
+  while (!nodeStack.empty()) 
+        {
+          BspNode *node = nodeStack.top();
+          nodeStack.pop();
+          
+          if (node->IsLeaf()) 
+                {
+                  ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
+                  
+                  if (!viewCell->Mailed()) 
                         {
-                                viewCell->Mail();
-                                viewCells.push_back(viewCell);
+                          viewCell->Mail();
+                          viewCells.push_back(viewCell);
                         }
                 }
                 else
                 {
-                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
-
-                        nodeStack.push(interior->GetFront());
-                        nodeStack.push(interior->GetBack());
+                  BspInterior *interior = dynamic_cast<BspInterior *>(node);
+                  
+                  nodeStack.push(interior->GetFront());
+                  nodeStack.push(interior->GetBack());
                 }
         }
@@ -2504 +2509 @@
 
 
+ViewCell *
+VspBspTree::GetViewCell(const Vector3 &point)
+{
+  if (mRoot == NULL)
+        return NULL;
+  
+  stack<BspNode *> nodeStack;
+  nodeStack.push(mRoot);
+  
+  ViewCell *viewcell = NULL;
+  
+  while (!nodeStack.empty())  {
+        BspNode *node = nodeStack.top();
+        nodeStack.pop();
+        
+        if (node->IsLeaf()) {
+          viewcell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
+          break;
+        } else {
+          
+          BspInterior *interior = dynamic_cast<BspInterior *>(node);
+                
+          // random decision
+          if (interior->GetPlane().Side(point) < 0)
+                nodeStack.push(interior->GetBack());
+          else
+                nodeStack.push(interior->GetFront());
+        }
+  }
+  
+  return viewcell;
+}
+
+
 int VspBspTree::RefineViewCells(const VssRayContainer &rays)
 {

trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.h

@@ -264 +264 @@
         BspNode *CollapseTree(BspNode *node);
 
+  ViewCell *
+  GetViewCell(const Vector3 &point);
+
         /** Constructs bsp rays for post processing and visualization.
         */

trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp

@@ -351 +351 @@
 {
 
-  mSceneGraph->CollectObjects(&mObjects);
 
   long startTime = GetTime();

trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.cpp

@@ -139 +139 @@
   point = GetOrigin() + t*dir;
   
-  return SqrDistance(point, center) < sqrRadius;
+  return ::SqrDistance(point, center) < sqrRadius;
 }
 
@@ -153 +153 @@
 {
   Vector3 dir = GetNormalizedDir();
+  return GetDirParam(axis, dir);
+}
+
+float
+VssRay::GetOpositeDirParametrization(const int axis) const
+{
+  Vector3 dir = -GetNormalizedDir();
   return GetDirParam(axis, dir);
 }
@@ -217 +224 @@
   
 }
+
+
+int
+VssRayContainer::SelectRays(const int number,
+                                                        VssRayContainer &selected) const
+{
+  float p = number/(float)size();
+  VssRayContainer::const_iterator it = begin();
+  
+  for (; it != end(); it++)
+        if (Random(1.0f) < p)
+          selected.push_back(*it);
+  
+  return selected.size();
+}
+
+
+int
+VssRayContainer::GetContributingRays(VssRayContainer &selected,
+                                                                         const int minPass
+                                                                         ) const
+{
+  VssRayContainer::const_iterator it = begin();
+  
+  for (; it != end(); it++) {
+        VssRay *ray = *it;
+        if (ray->mPass >= minPass && ray->mPvsContribution > 0)
+          selected.push_back(ray);
+  }
+  return selected.size();
+}
+
+

trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h

@@ -6 +6 @@
 #include "Vector3.h"
 #include "Ray.h"
+#include "Containers.h"
 
 class AxisAlignedBox3;
@@ -45 +46 @@
   Intersectable *mTerminationObject;
 
+  ViewCellContainer mViewCells;
+  
   ////////////////////////
   // members related to importance sampling
@@ -55 +58 @@
   // sum of relative ray contributions per object
   float mRelativePvsContribution;
+
+
+  // weighted contribution to the pvs (based on the pass the ray was casted at)
+  // computed by the prperocessor
+  float mWeightedPvsContribution;
   
   //////////////////////////////
@@ -71 +79 @@
         mRefCount(0),
     mFlags(0),
-        mPass(pass)
+        mPass(pass),
+        mViewCells(0),
+        mWeightedPvsContribution(0)
         //      mT(1.0f)
   {
@@ -82 +92 @@
         mMailbox(-1),
         mOriginObject(ray.sourceObject.mObject),
-        mPass(0)
+        mPass(0),
+        mViewCells(0)
         //      mT(1.0f)
-        {
-                if (ray.sourceObject.mObject)
-                        mOrigin = ray.Extrap(ray.sourceObject.mT);
-                else
-                        mOrigin = ray.GetLoc(); 
-
-                //Debug << "origin: " << mOrigin << endl;
-                if (!ray.intersections.empty())
-                {
-                        mTermination = ray.Extrap(ray.intersections[0].mT);
-                        mTerminationObject = ray.intersections[0].mObject;
-                }
-                else
-                {
-                        mTermination = ray.Extrap(1e6);//TODO: should be Limits::Infinity
-                        mTerminationObject = NULL;
-                }
-
-                Precompute();
-  }
-
+  {
+        if (ray.sourceObject.mObject)
+          mOrigin = ray.Extrap(ray.sourceObject.mT);
+        else
+          mOrigin = ray.GetLoc(); 
+        
+        //Debug << "origin: " << mOrigin << endl;
+        if (!ray.intersections.empty())
+          {
+                mTermination = ray.Extrap(ray.intersections[0].mT);
+                mTerminationObject = ray.intersections[0].mObject;
+          }
+        else
+          {
+                mTermination = ray.Extrap(1e6);//TODO: should be Limits::Infinity
+                mTerminationObject = NULL;
+          }
+        
+        Precompute();
+  }
+  
   void Precompute() {
     mFlags = 0;
@@ -158 +169 @@
         
   float GetDirParametrization(const int axis) const;
+  float GetOpositeDirParametrization(const int axis) const;
 
   static float VssRay::GetDirParam(const int axis, const Vector3 dir);
@@ -230 +242 @@
   }
 
+  friend bool GreaterWeightedPvsContribution(const VssRay * a,
+                                                                                         const VssRay *b) {
+        return a->mWeightedPvsContribution > b->mWeightedPvsContribution;
+  }
+
+  float SqrDistance(const Vector3 &point) const {
+        Vector3 diff = point - mOrigin;
+        float t = DotProd(diff, GetDir());
+        
+    if ( t <= 0.0f ) {
+          t = 0.0f;
+    } else {
+          if (t >= 1.0f) {
+                t = 1.0f;
+          } else {
+        t /= SqrMagnitude(GetDir());
+          }
+          diff -= t*GetDir();
+        }
+    return SqrMagnitude(diff);
+  }
+
 };
 
@@ -272 +306 @@
 {
   void PrintStatistics(ostream &s);
+  int SelectRays(const int number, VssRayContainer &selected) const;
+  int
+  GetContributingRays(VssRayContainer &selected,
+                                          const int minPass
+                                          ) const;
+  
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp

@@ -115 +115 @@
     else
       if ((*ri)->intersections.size()==0)
-        b = (*ri)->GetLoc() + length*(*ri)->GetDir();
+                b = (*ri)->GetLoc() + length*(*ri)->GetDir();
       else
-        b = (*ri)->Extrap((*ri)->intersections[0].mT);
+                b = (*ri)->Extrap((*ri)->intersections[0].mT);
     
     stream<<a.x<<" "<<a.y<<" "<<a.z<<" ,";

trunk/VUT/GtpVisibilityPreprocessor/src/default.env

@@ -15 +15 @@
 # filename ../data/atlanta/atlanta2.x3d
 # filename ../data/soda/soda.dat
- filename ../data/soda/soda5.dat
+filename ../data/soda/soda5.dat
 }
 
@@ -22 +22 @@
 #       type vss
         type rss
+        useGlRenderer true
 }
 
@@ -66 +67 @@
         initialSamples 500000
         vssSamples 10000000
-        vssSamplesPerPass 100000
+        vssSamplesPerPass 500000
         useImportanceSampling true
+
+        directionalSampling true
+        objectBasedSampling false
 
         Export {
@@ -77 +81 @@
 
         useViewcells true
-        updateSubdivision false
+        updateSubdivision true
 }
 
@@ -88 +92 @@
         minRays         30
         minSize         0.001
-        maxCostRatio    0.9
-        maxRayContribution 0.1
-        
+        maxCostRatio 1.0
+        maxRayContribution 0.5
+        maxRays         1000000
         maxTotalMemory  200
         maxStaticMemory 100
@@ -97 +101 @@
 #       splitType heuristic
         splitType hybrid
-        splitUseOnlyDrivingAxis false
+        splitUseOnlyDrivingAxis true
+        importanceBasedCost false
 
         interleaveDirSplits     true
@@ -157 +162 @@
         #type kdTree
         #type vspKdTree
-        #type bspTree
+        # type bspTree
         type vspBspTree
         
@@ -184 +189 @@
 }
 
-BspTree {
-        Construction {
-                samples 50000
-                epsilon 0.005
-        }
-
-        # random polygon       = 1
-        # axis aligned         = 2
-        # least splits         = 4
-        # balanced polygons    = 8
-        # balanced view cells  = 16
-        # largest polygon area = 32
-        # vertical axis        = 64
-        # blocked rays         = 128
-        # least ray splits     = 256
-        # balanced rays        = 512
-        # pvs                  = 1024
-
-        # least splits + balanced polygons
-        #splitPlaneStrategy 12
-        
-        #axis aligned + vertical axis
-        #splitPlaneStrategy 66
-        
-        # axis aligned + balanced view cells
-        # splitPlaneStrategy 18
-        
-        # largest polygon area
-        #splitPlaneStrategy 32
-        
-        # axus aligned + balanced polygons
-        #splitPlaneStrategy 72
-        
-        # axis aligned + blocked rays
-        #splitPlaneStrategy 130
-        
-        #splitPlaneStrategy 384
-        #splitPlaneStrategy 130
-        
-        splitPlaneStrategy 1024
-        
-        maxPolyCandidates 50
-        maxRayCandidates 50
-        
-        maxTests 10000
-        
-        # factors for evaluating split plane costs
-        Factor {
-                verticalSplits 1.0
-                largestPolyArea 1.0
-                blockedRays 1.0
-                leastRaySplits 1.0
-                balancedRays 1.0
-                pvs 1.0
-                leastSplits 1.0
-                balancedPolys 1.0
-                balancedViewCells 1.0
-        }
-        
-        Termination {
-                # parameters used for autopartition
-                minRays 200
-                minPolygons -1
-                maxDepth 40
-                minPvs 100
-                minArea 0.01
-                maxRayContribution 0.005
-                #maxAccRayLength 100
-                
-                # used for pvs criterium
-                ct_div_ci 0.0
-                
-                # axis aligned splits
-                AxisAligned {
-                        minPolys 5000
-                        minRays 500
-                        minObjects 10
-                        maxCostRatio 0.9
-                        ct_div_ci 0.5
-                }
-        }
-        
-        AxisAligned {
-                splitBorder 0.01
-        }
-        
-                
-        Visualization {
-        # x3d visualization of the split planes
-        exportSplits true
-        }
-}
 
 Simulation {
@@ -355 +268 @@
                         minRays 500
                         minObjects 10
-                        maxCostRatio 0.9
-                        ct_div_ci 0.5
+#                       maxCostRatio 0.9
+#                       ct_div_ci 0.5
                 }
         }
         
-        AxisAligned {
-                splitBorder 0.01
-        }
+#       AxisAligned {
+#               splitBorder 0.01
+#       }
         
         

trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp

@@ -13 +13 @@
 #include "ViewCell.h"
 #include "SceneGraph.h"
+#include "PreprocessorThread.h"
+
+#include <QApplication>
+#include <QtOpenGL>
+#include "GlRenderer.h"
 
 #define USE_EXE_PATH false
 
 
+
+
+
 int
-main(int argc, const char **argv)
+main(int argc, char **argv)
 {
   Debug.open("debug.log");
@@ -24 +32 @@
   environment->Parse(argc, argv, USE_EXE_PATH);
   MeshKdTree::ParseEnvironment();
-   
+  
   char buff[128];
   environment->GetStringValue("Preprocessor.type", buff);
@@ -48 +56 @@
                 }
   
+  preprocessor = p;
   
-        environment->GetStringValue("Scene.filename", buff);
-        string filename(buff);
-        
-        p->LoadScene(filename);
-        
-        p->BuildKdTree();
-        p->KdTreeStatistics(cout);
-        
-        // parse view cells related options
-        p->PrepareViewCells();
+  environment->GetStringValue("Scene.filename", buff);
+  string filename(buff);
+  p->LoadScene(filename);
+  
+  p->BuildKdTree();
+  p->KdTreeStatistics(cout);
+  
+  // parse view cells related options
+  p->PrepareViewCells();
+  
+
+  // create a preprocessor thread
+  PreprocessorThread *pt = new PreprocessorThread(p);
 
   //  p->mSceneGraph->Export("soda.x3d");
@@ -65 +77 @@
     p->Export(filename + "-kdtree.x3d", false, true, false);    
   }
+
+  if (p->mUseGlRenderer) {
+        QApplication a(argc, argv);
+        
+        if (!QGLFormat::hasOpenGL() || !QGLPbuffer::hasPbuffers()) {
+          QMessageBox::information(0, "OpenGL pbuffers", 
+                                                           "This system does not support OpenGL/pbuffers.",
+                                                           QMessageBox::Ok);
+          return -1;
+        }
+        
+        renderer = new GlRenderer(p->mSceneGraph, p->mViewCellsManager);
+        //  renderer->resize(640, 480);
+        renderer->resize(1024, 768);
+        renderer->show();
+        
+        
+        pt->start(QThread::LowPriority);
+        
+        return a.exec();
+  } else {
+        // just call the mail method -> will be executed in the main thread
+        pt->Main();
+  }
   
-   
-  if (1) {
-    p->ComputeVisibility();
-    p->ExportPreprocessedData("scene.vis");
-  }
-
-        Camera camera;
-  if (0) {
-    //camera.LookAtBox(p->mKdTree->GetBox());
-                camera.LookInBox(p->mKdTree->GetBox());
-                camera.SetPosition(camera.mPosition + Vector3(0,300,0));
-    camera.SnapImage("camera.jpg", p->mKdTree);
-        }
-        if (0) {
-    camera.LookInBox(p->mKdTree->GetBox());
-    camera.SetPosition(camera.mPosition - Vector3(0,100,0));
-    camera.SnapImage("camera2.jpg", p->mKdTree);
-  }
-
-        if (0) {
-    camera.SetPosition( p->mKdTree->GetBox().Center() - Vector3(0,-100,0) );
-                camera.SetDirection(Vector3(1, 0, 0));
-                camera.SnapImage("camera3.jpg", p->mKdTree);
-  }
-
   // clean up
-  DEL_PTR(p);
-  DEL_PTR(environment);
+  //  DEL_PTR(p);
+  //  DEL_PTR(environment);
 
   return 0;

trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro

@@ -16 +16 @@
 
 # debuc config
-CONFIG          += console warn_off thread release
+CONFIG          += console warn_off thread release qt
 
 # RELEASE CONFIG
@@ -23 +23 @@
 # DEPENDPATH    = ../../include
 
+QT += opengl
+
 win32:LIBS += xerces-c_2.lib devil.lib ilu.lib ilut.lib
 unix:LIBS += -lxerces-c -lIL -lILU -lILUT
 
 # Input
-HEADERS += Halton.h VssRay.h VssPreprocessor.h RssTree.h
+HEADERS += Halton.h VssRay.h VssPreprocessor.h RssTree.h GlRenderer.h \
+PreprocessorThread.h Preprocessor.h
 
 
@@ -39 +42 @@
 ViewCell.cpp ViewCellBsp.cpp Halton.cpp VssRay.cpp VssTree.cpp VssPreprocessor.cpp \
 RenderSimulator.cpp VspKdTree.cpp RayInfo.cpp RssTree.cpp RssPreprocessor.cpp \
-ViewCellsManager.cpp VspBspTree.cpp
+ViewCellsManager.cpp VspBspTree.cpp GlRenderer.cpp \
+PreprocessorThread.cpp Renderer.cpp