Changeset 382 for trunk

Timestamp:

11/05/05 20:03:25 (19 years ago)

Author:

bittner

Message:

vsspreprocessor update

Location:

trunk/VUT/GtpVisibilityPreprocessor/src

Files:

• Environment.cpp (modified) (1 diff)
• Intersectable.h (modified) (2 diffs)
• Makefile (modified) (7 diffs)
• Mesh.cpp (modified) (1 diff)
• VssPreprocessor.cpp (modified) (11 diffs)
• VssPreprocessor.h (modified) (1 diff)
• VssRay.cpp (modified) (1 diff)
• VssRay.h (modified) (2 diffs)
• VssTree.cpp (modified) (24 diffs)
• VssTree.h (modified) (14 diffs)
• default.env (modified) (3 diffs)
• preprocessor.pro (modified) (1 diff)

trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp

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

trunk/VUT/GtpVisibilityPreprocessor/src/Intersectable.h

@@ -11 +11 @@
         // both mailId and mailbox should be unique for each thread!!!
         static int sMailId;
+        static int sReservedMailboxes;
         int mMailbox;
         
@@ -26 +27 @@
         int GetId() { return mId; }
         
+
+        static void NewMail(const int reserve = 1) {
+                sMailId += sReservedMailboxes;
+                sReservedMailboxes = reserve;
+        }
+        
         void Mail() { mMailbox = sMailId; }
-  static void NewMail() { sMailId++; }
   bool Mailed() const { return mMailbox == sMailId; }
-  int IncMail() { return ++mMailbox - sMailId; }
 
+        void Mail(const int mailbox) { mMailbox = sMailId + mailbox; }
+        bool Mailed(const int mailbox) const { return mMailbox == sMailId + mailbox; }
+
+        int IncMail() { return ++mMailbox - sMailId; }
+
+
+        
   virtual AxisAlignedBox3 GetBox() = 0;
   

trunk/VUT/GtpVisibilityPreprocessor/src/Makefile

@@ -1 +1 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (1.07a) (Qt 3.3.2) on: Thu Nov 03 11:34:08 2005
+# Generated by qmake (1.07a) (Qt 3.3.2) on: Sat Nov 05 19:00:35 2005
 # Project:  preprocessor.pro
 # Template: app
@@ -415 +415 @@
 main.obj: main.cpp  \
                 SamplingPreprocessor.h \
+                VssPreprocessor.h \
                 ExactPreprocessor.h \
                 Parser.h \
@@ -441 +442 @@
                 common.h \
                 Ray.h \
+                VssRay.h \
                 
 
@@ -639 +641 @@
                 Environment.h \
                 VssRay.h \
+                Intersectable.h \
                 AxisAlignedBox3.h \
                 Statistics.h \
@@ -646 +649 @@
                 Matrix4x4.h \
                 Plane3.h \
+                Pvs.h \
                 
 
@@ -657 +661 @@
                 Polygon3.h \
                 ViewCell.h \
+                VssRay.h \
+                VssTree.h \
                 Containers.h \
                 AxisAlignedBox3.h \
@@ -672 +678 @@
                 Material.h \
                 Exporter.h \
+                Statistics.h \
                 
 

trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp

@@ -5 +5 @@
 
 int Intersectable::sMailId = 21843194198;
+int Intersectable::sReservedMailboxes = 1;
 
 void

trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp

@@ -8 +8 @@
 #include "ViewCell.h"
 #include "VssRay.h"
+#include "VssTree.h"
 
 VssPreprocessor::VssPreprocessor():
@@ -48 +49 @@
         Intersectable *objectA, *objectB;
         Vector3 pointA, pointB;
-
+        float bsize = Magnitude(box.Size());
         if (mKdTree->CastRay(ray)) {
                 objectA = ray.intersections[0].mObject;
@@ -57 +58 @@
                 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);
+                
         }
         
@@ -69 +75 @@
                 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);
         }
@@ -86 +97 @@
 
 Vector3
-VssPreprocessor::GetViewpoint()
-{
-        AxisAlignedBox3 box = mKdTree->GetBox();
-
+VssPreprocessor::GetViewpoint(AxisAlignedBox3 *viewSpaceBox)
+{
+        AxisAlignedBox3 box;
+        
+        if (viewSpaceBox)
+                box =*viewSpaceBox;
+        else
+                box = mKdTree->GetBox();
+        
         // shrink the box in the y direction
-        Vector3 diff(0, -box.Size().y*0.4f, 0);
-        box.Enlarge(diff);
-        
         return box.GetRandomPoint();
 }
@@ -100 +113 @@
 VssPreprocessor::GetDirection(const Vector3 &viewpoint)
 {
-        int i = RandomValue(0, mObjects.size());
+        int i = RandomValue(0, mObjects.size()-1);
         Intersectable *object = mObjects[i];
         Vector3 point, normal;
@@ -114 +127 @@
   
   mSceneGraph->CollectObjects(&mObjects);
-        
-
+        cout<<"#NUM_OBJECTS (Total numner of objects)\n"<<mObjects.size()<<endl;
+        
   long startTime = GetTime();
   
@@ -121 +134 @@
   int totalSamples = 0;
 
-
+        AxisAlignedBox3 *viewSpaceBox = NULL;
+
+        AxisAlignedBox3 box = mKdTree->GetBox();
+        box.Enlarge(box.Size()*-Vector3(0.45, 0.45, 0.45));
+
+        bool useViewSpaceBox = false;
+        if (useViewSpaceBox)
+                viewSpaceBox = &box;
+        
   while (totalSamples < mTotalSamples) {
                 int passContributingSamples = 0;
@@ -132 +153 @@
                 
                 for (int k=0; k < mSamplesPerPass; k++) {
-
-                        Vector3 viewpoint = GetViewpoint();
+                        
+                        Vector3 viewpoint = GetViewpoint(viewSpaceBox);
                         Vector3 direction = GetDirection(viewpoint);
-
+                        
                         VssRay *vssRay = CastRay(viewpoint, direction);
                         
@@ -153 +174 @@
     
                 mPass++;
-
+                
                 int pvsSize = 0;
                 float avgRayContrib = (passContributingSamples > 0) ? 
@@ -176 +197 @@
         
         cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;
-        cout << "#totalRayStackSize=" << mVssRays.size() << endl;
-
+        cout << "#totalRayStackSize=" << mVssRays.size() << endl <<flush;
+
+        VssTree *vssTree = new VssTree;
+
+        vssTree->Construct(mVssRays, viewSpaceBox);
+        
   return true;
 }

trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.h

@@ -28 +28 @@
 
         Vector3
-        GetViewpoint();
+        GetViewpoint(AxisAlignedBox3 *viewSpaceBox);
 
         Vector3

trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.cpp

@@ -141 +141 @@
   return SqrDistance(point, center) < sqrRadius;
 }
+
+
+float
+VssRay::GetDirParametrization(const int axis) const
+{
+        float angle;
+        Vector3 dir = GetNormalizedDir();
+        if (axis == 0) {
+                angle = atan2(dir.x, dir.z);
+        } else {
+                angle = asin(dir.y);
+        }
+        return angle;
+}

trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h

@@ -22 +22 @@
         
         // side of the ray - used for the ray classification
-  char mSide;
+        //  char mSide;
         
         // computed t
@@ -98 +98 @@
   Vector3 GetNormalizedDir() const { return (mTermination - mOrigin)*mInvSize; }
 
+        float GetDirParametrization(const int axis) const;
+        
         float GetInvSize() const { return  mInvSize; }
   float GetOrigin(const int axis) const { return mOrigin[axis]; }

trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp

@@ -22 +22 @@
 #include "Environment.h"
 #include "VssRay.h"
+#include "Intersectable.h"
 
 
@@ -32 +33 @@
 VssTreeLeaf::mailID = 0;
 
+inline void
+AddObject2Pvs(Intersectable *object,
+                                                        const int side,
+                                                        int &pvsBack,
+                                                        int &pvsFront)
+{
+
+        if (!object)
+                return;
+                
+        if (side <= 0) {
+                if (!object->Mailed() && !object->Mailed(2)) {
+                        pvsBack++;
+                        if (object->Mailed(1))
+                                object->Mail(2);
+                        else
+                                object->Mail();
+                }
+        }
+        
+        if (side >= 0) {
+                if (!object->Mailed(1) && !object->Mailed(2)) {
+                        pvsFront++;
+                        if (object->Mailed())
+                                object->Mail(2);
+                        else
+                                object->Mail(1);
+                }
+        }
+}
+
 
 // Constructor
@@ -37 +69 @@
 {
   environment->GetIntValue("VssTree.maxDepth", termMaxDepth);
-
-        environment->GetIntValue("VssTree.minCost", termMinCost);
+        environment->GetIntValue("VssTree.minPvs", termMinPvs);
+        environment->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution);
+  environment->GetFloatValue("VssTree.maxCostRatio", termMaxCostRatio);
 
   environment->GetFloatValue("VssTree.minSize", termMinSize);
   termMinSize = sqr(termMinSize);
-
+        
   environment->GetFloatValue("VssTree.refDirBoxMaxSize", refDirBoxMaxSize);
   refDirBoxMaxSize = sqr(refDirBoxMaxSize);
@@ -48 +81 @@
   environment->GetFloatValue("VssTree.epsilon", epsilon);
   environment->GetFloatValue("VssTree.ct_div_ci", ct_div_ci);
-
+        
   environment->GetFloatValue("VssTree.maxTotalMemory", maxTotalMemory);
   environment->GetFloatValue("VssTree.maxStaticMemory", maxStaticMemory);
   
-  environment->GetFloatValue("VssTree.maxCostRatio", maxCostRatio);
   
 
@@ -77 +109 @@
     splitType = ESplitRegular;
   else
-    if (name.compare("volume") == 0)
+    if (name.compare("heuristics") == 0)
       splitType = ESplitVolume;
-    else
-      if (name.compare("queries") == 0)
-                                splitType = ESplitQueries;
-      else {
-                                cerr<<"Invalid VssTree split type "<<name<<endl;
-                                exit(1);
-      }
+                else {
+                        cerr<<"Invalid VssTree split type "<<name<<endl;
+                        exit(1);
+                }
         
   environment->GetBoolValue("VssTree.randomize", randomize);
@@ -91 +120 @@
   root = NULL;
   
-  splitCandidates = new vector<SSortableEntry>;
+  splitCandidates = new vector<SortableEntry>;
 }
 
@@ -109 +138 @@
   app << "===== VssTree statistics ===============\n";
 
-  app << "#N_RAYS Number of rays )\n"
+  app << "#N_RAYS ( Number of rays )\n"
       << rays <<endl;
-  app << "#N_DOMAINS  ( Number of query domains )\n"
-      << queryDomains <<endl;
+
+        app << "#N_INITPVS ( Initial PVS size )\n"
+      << initialPvsSize <<endl;
   
   app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
@@ -135 +165 @@
     maxRayRefs << "\n";
 
-  app << "#N_NONEMPTYRAYREFS  ( Number of rayRefs in nonEmpty leaves / non empty leaf )\n" <<
-    rayRefsNonZeroQuery/(double)(Leaves() - zeroQueryNodes) << "\n";
-
-  app << "#N_LEAFDOMAINREFS  ( Number of query domain Refs / leaf )\n" <<
-    queryDomainRefs/(double)Leaves() << "\n";
 
   //  app << setprecision(4);
-
-  app << "#N_PEMPTYLEAVES  ( Percentage of leaves with zero query domains )\n"<<
-    zeroQueryNodes*100/(double)Leaves()<<endl;
 
   app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
@@ -151 +173 @@
   app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
     minCostNodes*100/(double)Leaves()<<endl;
+        
+        app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
+    minSizeNodes*100/(double)Leaves()<<endl;
+
+        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
+    maxRayContribNodes*100/(double)Leaves()<<endl;
 
   app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
@@ -173 +201 @@
 VssTree::Construct(
                                                                          VssRayContainer &rays,
-                                                                         const bool onlyStaticPart
+                                                                         AxisAlignedBox3 *forcedBoundingBox
                                                                          )
 {
   stat.Start();
   
-  if (onlyStaticPart)
-    maxMemory = maxStaticMemory;
-  else
-    maxMemory = maxTotalMemory;
-
+        maxMemory = maxStaticMemory;
 
   if (root)
@@ -194 +218 @@
   
   bbox.Initialize();
-
   dirBBox.Initialize();
   
@@ -201 +224 @@
       ri++) {
     leaf->AddRay(VssTreeNode::RayInfo(*ri));
+
     bbox.Include((*ri)->GetOrigin());
     bbox.Include((*ri)->GetTermination());
     
-    dirBBox.Include( (*ri)->GetNormalizedDir() );
-  }
-  
-  
+                
+                dirBBox.Include(Vector3(
+                                                                                                                (*ri)->GetDirParametrization(0),
+                                                                                                                (*ri)->GetDirParametrization(1),
+                                                                                                                0
+                                                                                                                )
+                                                                                );
+        }
+        
+        
+        if (forcedBoundingBox)
+                bbox = *forcedBoundingBox;
+
+        cout<<"Bbox = "<<bbox<<endl;
+        cout<<"Dirr Bbox = "<<dirBBox<<endl;
+
   stat.rays = leaf->rays.size();
-  
+        EvalLeafPvs(leaf);
+  stat.initialPvsSize = leaf->mPvsSize;
   // Subdivide();
   root = Subdivide(TraversalData(leaf, bbox, 0));
@@ -216 +253 @@
     // force realease of this vector
     delete splitCandidates;
-    splitCandidates = new vector<SSortableEntry>;
+    splitCandidates = new vector<SortableEntry>;
   }
   
   stat.Stop();
 
-  cout<<"##################################"<<endl;
-  
+        stat.Print(cout);
+        cout<<"#Total memory="<<GetMemUsage()<<endl;
+
+}
+
+void
+VssTree::EvalLeafPvs(VssTreeLeaf *leaf)
+{
+        Intersectable::NewMail();
+        int pvsSize = 0;
+        for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                        ri != leaf->rays.end();
+                        ri++)
+                if ((*ri).mRay->IsActive()) {
+                        Intersectable *object = (*ri).mRay->mOriginObject;
+                        if (object && !object->Mailed()) {
+                                pvsSize++;
+                                object->Mail();
+                        }
+                        object = (*ri).mRay->mTerminationObject;
+                        if (object && !object->Mailed()) {
+                                pvsSize++;
+                                object->Mail();
+                        }
+                }
+        leaf->mPvsSize = pvsSize;
 }
 
@@ -283 +344 @@
                                                                                  VssTreeLeaf *leaf,
                                                                                  const AxisAlignedBox3 &box,
-                                                                                 float &position
+                                                                                 float &position,
+                                                                                 int &raysBack,
+                                                                                 int &raysFront,
+                                                                                 int &pvsBack,
+                                                                                 int &pvsFront
                                                                                  )
 {
 
   int minDirDepth = 6;
-
-  Vector3 size = box.Size();
-  
-  float pDirSubdivision = 0.0f;
-  if (leaf->depth >= minDirDepth && RandomValue(0, 1) < pDirSubdivision) {
-    AxisAlignedBox3 dirbox = GetDirBBox(leaf);
-    Vector3 size = dirbox.Size();
-    int axis = size.DrivingAxis();
-    position = (dirbox.Min()[axis] + dirbox.Max()[axis])*0.5f;
-    return axis + 3;
-  }
-  
-  int axis = size.DrivingAxis();
-
-  
+  int axis;
+        float costRatio;
+        
   if (splitType == ESplitRegular) {
-                position = (box.Min()[axis] + box.Max()[axis])*0.5f;
-  }
-  
-  return axis;
+                int sRaysBack, sRaysFront;
+                int dRaysBack, dRaysFront;
+                int sPvsBack, sPvsFront;
+                int dPvsBack, dPvsFront;
+                
+                
+                int sAxis = box.Size().DrivingAxis();
+                float sPosition = (box.Min()[sAxis] + box.Max()[sAxis])*0.5f;
+                float sCostRatio = EvalCostRatio(leaf,
+                                                                                                                                                 sAxis,
+                                                                                                                                                 sPosition,
+                                                                                                                                                 sRaysBack,
+                                                                                                                                                 sRaysFront,
+                                                                                                                                                 sPvsBack,
+                                                                                                                                                 sPvsFront
+                                                                                                                                                 );
+                //              cout<<"srays back="<<sRaysBack<<" rays front="<<sRaysFront<<" pvs back="<<sPvsBack<<
+                //                      " pvs front="<<sPvsFront<<endl;
+
+                AxisAlignedBox3 dirBox = GetDirBBox(leaf);
+                int dAxis = dirBox.Size().DrivingAxis();
+                float dPosition = (dirBox.Min()[dAxis] + dirBox.Max()[dAxis])*0.5f;
+                float dCostRatio = EvalCostRatio(leaf,
+                                                                                                                                                 dAxis+3,
+                                                                                                                                                 dPosition,
+                                                                                                                                                 dRaysBack,
+                                                                                                                                                 dRaysFront,
+                                                                                                                                                 dPvsBack,
+                                                                                                                                                 dPvsFront
+                                                                                                                                                 );
+                
+                //              cout<<"drays back="<<dRaysBack<<" rays front="<<dRaysFront<<" pvs back="<<dPvsBack<<
+                //                      " pvs front="<<dPvsFront<<endl;
+
+                if (sCostRatio < dCostRatio) {
+                        costRatio = sCostRatio;
+                        axis = sAxis;
+                        position = sPosition;
+                        raysBack = sRaysBack;
+                        raysFront = sRaysFront;
+                        pvsBack = sPvsBack;
+                        pvsFront = sPvsFront;
+                        
+                } else {
+                        costRatio = dCostRatio;
+                        axis = dAxis+3;
+                        position = dPosition;
+                        raysBack = dRaysBack;
+                        raysFront = dRaysFront;
+                        pvsBack = dPvsBack;
+                        pvsFront = dPvsFront;
+                }
+                
+        } else {
+    if (splitType == ESplitVolume)
+      costRatio = BestCostRatio(leaf,
+                                                                                                                                axis,
+                                                                                                                                position,
+                                                                                                                                raysBack,
+                                                                                                                                raysFront);
+                else {
+                        cerr<<"VssTree: Unknown split heuristics\n";
+                        exit(1);
+                }
+  }
+
+        if (costRatio > termMaxCostRatio) {
+                cout<<"Too big cost ratio "<<costRatio<<endl;
+                return -1;
+        }
+
+#if 0   
+        cout<<
+                "pvs="<<leaf->mPvsSize<<
+                " rays="<<leaf->rays.size()<<
+                " rc="<<leaf->GetAvgRayContribution()<<
+                " axis="<<axis<<endl;
+#endif
+        
+        return axis;
+}
+
+
+                                                        
+
+float
+VssTree::EvalCostRatio(
+                                                                                         VssTreeLeaf *leaf,
+                                                                                         const int axis,
+                                                                                         const float position,
+                                                                                         int &raysBack,
+                                                                                         int &raysFront,
+                                                                                         int &pvsBack,
+                                                                                         int &pvsFront
+                                                                                         )
+{
+        raysBack = 0;
+        raysFront = 0;
+        pvsFront = 0;
+        pvsBack = 0;
+
+        float newCost;
+
+        Intersectable::NewMail(3);
+        
+        // eval pvs size
+        int pvsSize = leaf->mPvsSize;
+
+        Intersectable::NewMail(3);
+
+        if (axis <= VssTreeNode::SPLIT_Z) {
+    // this is the main ray classification loop!
+    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                                ri != leaf->rays.end();
+                                ri++)
+      if ((*ri).mRay->IsActive()) {
+                                
+                                // determine the side of this ray with respect to the plane
+                                int side = (*ri).ComputeRayIntersection(axis, position, (*ri).mRay->mT);
+                                
+                                //                              (*ri).mRay->mSide = side;
+                                
+                                if (side <= 0)
+                                        raysBack++;
+                                
+                                if (side >= 0)
+                                        raysFront++;
+                                
+                                AddObject2Pvs((*ri).mRay->mOriginObject, side, pvsBack, pvsFront);
+                                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+
+      }
+                
+                AxisAlignedBox3 box = GetBBox(leaf);
+                
+                float minBox = box.Min(axis);
+                float maxBox = box.Max(axis);
+                float sizeBox = maxBox - minBox;
+                
+                //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+                float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
+
+                newCost = ct_div_ci + sum/sizeBox;
+                
+  } else {
+                
+                // directional split
+    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                                ri != leaf->rays.end();
+                                ri++)
+      if ((*ri).mRay->IsActive()) {
+                                
+                                // determine the side of this ray with respect to the plane
+                                int side;
+                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
+                                        side = 1;
+                                else
+                                        side = -1;
+
+                                if (side <= 0)
+                                        raysBack++;
+                                
+                                if (side >= 0)
+                                        raysFront++;
+
+                                //                              (*ri).mRay->mSide = side;
+
+                                AddObject2Pvs((*ri).mRay->mOriginObject, side, pvsBack, pvsFront);
+                                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
+
+      }
+                
+                AxisAlignedBox3 box = GetDirBBox(leaf);
+                
+                float minBox = box.Min(axis-3);
+                float maxBox = box.Max(axis-3);
+                float sizeBox = maxBox - minBox;
+                
+                //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
+                float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
+                //              float sum = pvsBack + pvsFront;
+                newCost = ct_div_ci + sum/sizeBox;
+  }
+
+        //      cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
+        //  float oldCost = leaf->rays.size();
+        float oldCost = pvsSize;
+  float ratio = newCost/oldCost;
+        //      cout<<"ratio="<<ratio<<endl;
+        return ratio;
+}
+
+float
+VssTree::BestCostRatio(
+                                                                                         VssTreeLeaf *node,
+                                                                                         int &axis,
+                                                                                         float &position,
+                                                                                         int &raysBack,
+                                                                                         int &raysFront
+                                                                                         )
+{
+        AxisAlignedBox3 box = GetBBox(node);
+        AxisAlignedBox3 dirBox = GetDirBBox(node);
+
+        axis = box.Size().DrivingAxis();
+        
+        SortSplitCandidates(node, axis);
+  
+  // go through the lists, count the number of objects left and right
+  // and evaluate the following cost funcion:
+  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
+
+  int rl=0, rr=node->rays.size();
+        
+  float minBox = box.Min(axis);
+  float maxBox = box.Max(axis);
+  float sizeBox = maxBox - minBox;
+  
+  float minBand = minBox + 0.1*(maxBox - minBox);
+  float maxBand = minBox + 0.9*(maxBox - minBox);
+  
+  float sum = rr*sizeBox;
+  float minSum = 1e20;
+  
+  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
+      ci < splitCandidates->end();
+      ci++) {
+    
+    switch ((*ci).type) {
+    case SortableEntry::ERayMin:
+      rl++;
+      break;
+    case SortableEntry::ERayMax:
+      rr--;
+      break;
+    }
+                
+    if ((*ci).value > minBand && (*ci).value < maxBand) {
+      
+      sum = rl*((*ci).value - minBox) + rr*(maxBox - (*ci).value);
+      
+      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
+      //      cout<<"cost= "<<sum<<endl;
+      
+      if (sum < minSum) {
+                                minSum = sum;
+                                position = (*ci).value;
+                                
+                                raysBack = rl;
+                                raysFront = rr;
+                                
+      }
+    }
+  }
+  
+  float oldCost = node->rays.size();
+  float newCost = ct_div_ci + minSum/sizeBox;
+  float ratio = newCost/oldCost;
+  
+  //  cout<<"===================="<<endl;
+  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
+  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
+
+        return ratio;
+}
+
+void
+VssTree::SortSplitCandidates(
+                                                                                                                 VssTreeLeaf *node,
+                                                                                                                 const int axis
+                                                                                                                 )
+{
+  
+  splitCandidates->clear();
+  
+  int requestedSize = 2*(node->rays.size());
+  // creates a sorted split candidates array
+  if (splitCandidates->capacity() > 500000 &&
+      requestedSize < (int)(splitCandidates->capacity()/10) ) {
+    
+    delete splitCandidates;
+    splitCandidates = new vector<SortableEntry>;
+  }
+  
+  splitCandidates->reserve(requestedSize);
+
+  // insert all queries 
+  for(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
+      ri < node->rays.end();
+      ri++) {
+    bool positive = (*ri).mRay->HasPosDir(axis);
+    splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
+                                                                                                                                                                                 SortableEntry::ERayMax,
+                                                                                                                                                                                 (*ri).ExtrapOrigin(axis),
+                                                                                                                                                                                 (void *)&*ri)
+                                                                                                                         );
+                
+    splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
+                                                                                                                                                                                 SortableEntry::ERayMin,
+                                                                                                                                                                                 (*ri).ExtrapTermination(axis),
+                                                                                                                                                                                 (void *)&*ri)
+                                                                                                                         );
+  }
+        
+  stable_sort(splitCandidates->begin(), splitCandidates->end());
 }
 
@@ -318 +672 @@
   VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
 
-  if (data.depth > termMaxDepth)
+  if (data.depth >= termMaxDepth)
     stat.maxDepthNodes++;
   
-  if ( (int)(leaf->rays.size()) < termMinCost)
-    stat.minCostNodes++;
-  
-
-  if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
+        //  if ( (int)(leaf->rays.size()) < termMinCost)
+        //    stat.minCostNodes++;
+        if ( leaf->mPvsSize < termMinPvs)
+                stat.minCostNodes++;
+  
+        if (leaf->GetAvgRayContribution() > termMaxRayContribution )
+                stat.maxRayContribNodes++;
+        
+        if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
+                stat.minSizeNodes++;
+        }
+
+        if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
     stat.maxRayRefs = leaf->rays.size();
-  
+
 }
 
@@ -341 +703 @@
 {
   
-  if ( ((int)leaf->rays.size() < termMinCost) ||
-       (leaf->depth >= termMaxDepth) ) {
-    return leaf;
-  }
-  
-  float position;
-  
+  if ( (leaf->mPvsSize < termMinPvs) ||
+                         (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
+       (leaf->depth >= termMaxDepth) ||
+                         SqrMagnitude(box.Size()) <= termMinSize ) {
+
+#if 0
+                if (leaf->depth >= termMaxDepth) {
+                        cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
+                        cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
+                }
+#endif
+                
+                return leaf;
+        }
+        
+        float position;
+        
+        // first count ray sides
+  int raysBack;
+  int raysFront;
+        int pvsBack;
+        int pvsFront;
+        
   // select subdivision axis
-  int axis = SelectPlane( leaf, box, position);
+  int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
+
+        //      cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
+        //              pvsFront<<endl;
 
   if (axis == -1) {
@@ -364 +745 @@
   node->position = position;
   node->bbox = box;
-
   node->dirBBox = GetDirBBox(leaf);
   
@@ -370 +750 @@
   frontBBox = box;
 
-  // first count ray sides
-  int raysBack = 0;
-  int raysFront = 0;
-  
-  if (axis <= VssTreeNode::SPLIT_Z) {
-    backBBox.SetMax(axis, position);
-    frontBBox.SetMin(axis, position);
-    
-    // this is the main ray classification loop!
-    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                                ri != leaf->rays.end();
-                                ri++)
-      if ((*ri).mRay->IsActive()) {
-                                
-                                // determine the side of this ray with respect to the plane
-                                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
-                                
-                                (*ri).mRay->mSide = side;
-                                
-                                if (side <= 0)
-                                        raysBack++;
-                                
-                                if (side >= 0)
-                                        raysFront++;
-                                
-#if 0
-        cout<<"-------------------"<<endl;
-        cout<<"axis = "<<(int)node->axis<<"\t position="<<node->position<<endl;
-        cout<<"origin="<<(*ri).ray->GetOrigin()<<"\t  term="<<(*ri).ray->GetTermination()<<endl;
-        cout<<"side = "<<side<<"\t t="<<t<<endl;
-        cout<<"-------------------"<<endl;
-#endif
-      }
-  } else {
-    // directional split
-    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
-                                ri != leaf->rays.end();
-                                ri++)
-      if ((*ri).mRay->IsActive()) {
-                                
-                                // determine the side of this ray with respect to the plane
-                                int side;
-                                if ((*ri).mRay->GetNormalizedDir(axis - 3) > position)
-                                        side = 1;
-                                else
-                                        side = -1;
-                                
-                                
-                                (*ri).mRay->mSide = side;
-                                
-                                if (side <= 0)
-                                        raysBack++;
-                                
-                                if (side >= 0)
-                                        raysFront++;
-      }
-  }
-  
   VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
+        back->mPvsSize = pvsBack;
   VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
+        front->mPvsSize = pvsFront;
 
   // replace a link from node's parent
@@ -436 +760 @@
   // and setup child links
   node->SetupChildLinks(back, front);
-
-  
+        
   if (axis <= VssTreeNode::SPLIT_Z) {
-    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
+                backBBox.SetMax(axis, position);
+    frontBBox.SetMin(axis, position);
+                
+                for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
                                 ri != leaf->rays.end();
                                 ri++) {
@@ -446 +772 @@
                                 // first unref ray from the former leaf
                                 (*ri).mRay->Unref();
-                                
-                                if ((*ri).mRay->mSide == 0) {
+
+                                // determine the side of this ray with respect to the plane
+                                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
+
+                                if (side == 0) {
                                         if ((*ri).mRay->HasPosDir(axis)) {
                                                 back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
@@ -465 +794 @@
                                         }
                                 } else
-                                        if ((*ri).mRay->mSide == 1) 
+                                        if (side == 1) 
                                                 front->AddRay(*ri);
                                         else
@@ -475 +804 @@
     // rays front/back
     
-#if DEBUG_DIR_SPLIT
-    cout<<"dir split, depth="<<(int)leaf->depth<<" front= "<<raysFront<<" back="<<raysBack<<endl;
-#endif
     
     for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
@@ -485 +811 @@
                                 // first unref ray from the former leaf
                                 (*ri).mRay->Unref();
-                                
-                                if ((*ri).mRay->mSide == 1)
+
+                                int side;
+                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
+                                        side = 1;
+                                else
+                                        side = -1;
+                                
+                                if (side == 1)
                                         front->AddRay(*ri);
                                 else
@@ -573 +905 @@
   // check if we should perform a dynamic subdivision of the leaf
   if (
-      leaf->rays.size() > (unsigned)termMinCost &&
+                        //      leaf->rays.size() > (unsigned)termMinCost &&
+                        (leaf->mPvsSize >= termMinPvs) &&
       SqrMagnitude(leafBBox.Size()) > sizeThreshold) {
     

trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.h

@@ -43 +43 @@
   // totals number of rays
   int rays;
-  // total number of query domains
+        // initial size of the pvs
+  int initialPvsSize;
+        // total number of query domains
   int queryDomains;
   // total number of ray references
   int rayRefs;
-  // refs in non empty leafs
-  int rayRefsNonZeroQuery;
-  // total number of query references
-  int queryDomainRefs;
-  // nodes with zero queries
-  int zeroQueryNodes;
-  // max depth nodes
+
+        // max depth nodes
   int maxDepthNodes;
   // max depth nodes
   int minCostNodes;
+        // max ray contribution nodes
+  int maxRayContribNodes;
+        // max depth nodes
+  int minSizeNodes;
+        
   // max number of rays per node
   int maxRayRefs;
@@ -78 +80 @@
       splits[i] = 0;
     rays = queryDomains = 0;
-    rayRefs = rayRefsNonZeroQuery = queryDomainRefs = 0;
-    zeroQueryNodes = 0;
+    rayRefs = 0;
     maxDepthNodes = 0;
     minCostNodes = 0;
     maxRayRefs = 0;
     addedRayRefs = removedRayRefs = 0;
+                initialPvsSize = 0;
+                maxRayContribNodes = 0;
+                minSizeNodes = 0;
   }
 
@@ -101 +105 @@
 // For sorting rays
 // --------------------------------------------------------------
-struct  SSortableEntry
+struct  SortableEntry
 {
   enum EType {
@@ -112 +116 @@
   void *data;
   
-  SSortableEntry() {}
-  SSortableEntry(const int t, const float v, void *d):type(t),
-                                                                                                                                                                                                                        value(v),
-                                                      data(d) {}
-        
-  friend bool operator<(const SSortableEntry &a, const SSortableEntry &b) {
+  SortableEntry() {}
+  SortableEntry(const int t, const float v, void *d):type(t),
+                                                                                                                                                                                                                 value(v),
+                                                                                                                                                                                                                 data(d) {}
+        
+  friend bool operator<(const SortableEntry &a, const SortableEntry &b) {
     return a.value < b.value;
   }
@@ -215 +219 @@
         void SetMaxT (const float t) { mMaxT = t; }
 #endif 
+
+
+  int ComputeRayIntersection(const float axis,
+                                                                                                                 const float position,
+                                                                                                                 float &t
+                                                                                                                 ) const {
+                
+                // intersect the ray with the plane
+    float denom = mRay->GetDir(axis);
+    
+    if (fabs(denom) < 1e-20)
+      //if (denom == 0.0f)
+      return (mRay->GetOrigin(axis) > position) ? 1 : -1;
+    
+    t = (position - mRay->GetOrigin(axis))/denom;
+
+    if (t < GetMinT())
+      return (denom > 0) ? 1 : -1;
+
+    if (t > GetMaxT())
+      return (denom > 0) ? -1 : 1;
+
+                return 0;
+        }
+
+
 };
 
@@ -329 +359 @@
 
   
+        
   int ComputeRayIntersection(const RayInfo &rayData,
                                                                                                                  float &t
                                                                                                                  ) {
-                
-    // intersect the ray with the plane
-    float denom = rayData.mRay->GetDir(axis);
-    
-    if (fabs(denom) < 1e-20)
-      //if (denom == 0.0f)
-      return (rayData.mRay->GetOrigin(axis) > position) ? 1 : -1;
-    
-    t = (position - rayData.mRay->GetOrigin(axis))/denom;
-
-    if (t < rayData.GetMinT())
-      return (denom > 0) ? 1 : -1;
-
-    if (t > rayData.GetMaxT())
-      return (denom > 0) ? -1 : 1;
-    
-    return 0;
+                return rayData.ComputeRayIntersection(axis, position, t);
   }
 
@@ -365 +380 @@
   
   RayInfoContainer rays;
-
-  
-  VssTreeLeaf(VssTreeInterior *p, const int n):VssTreeNode(p), rays() {
-    rays.reserve(n);
+        int mPvsSize;
+  
+  VssTreeLeaf(VssTreeInterior *p,
+                                                        const int nRays
+                                                        ):VssTreeNode(p), rays(), mPvsSize(0) {
+    rays.reserve(nRays);
   }
   
@@ -392 +409 @@
   }
 
+        float GetAvgRayContribution() const {
+                return mPvsSize/(float)rays.size();
+        }
 };
 
@@ -490 +510 @@
   // max depth of the tree
   int termMaxDepth;
-  // minimal cost of the node to still get subdivided
-  int termMinCost;
   // minimal ratio of the volume of the cell and the query volume
   float termMinSize;
 
-  // minimal cost imporvement to subdivide a node
-  float maxCostRatio;
-  
+        // minimal pvs per node to still get subdivided
+  int termMinPvs;
+
+  // maximal cost ration to subdivide a node
+  float termMaxCostRatio;
+        
+        // maximal contribution per ray to subdivide the node
+        float termMaxRayContribution;
+
+        
   // randomized construction
   bool randomize;
 
   // type of the splitting to use fo rthe tree construction
-  enum {ESplitRegular, ESplitVolume, ESplitQueries };
+  enum {ESplitRegular, ESplitVolume };
   int splitType;
         
@@ -528 +553 @@
   
   // reusable array of split candidates
-  vector<SSortableEntry> *splitCandidates;
+  vector<SortableEntry> *splitCandidates;
   /////////////////////////////
 
@@ -540 +565 @@
   Construct(
                                                 VssRayContainer &rays,
-                                                const bool onlyStaticPart
+                                                AxisAlignedBox3 *forcedBoundingBox = NULL
                                                 );
         
@@ -565 +590 @@
   int
   SelectPlane(VssTreeLeaf *leaf,
-              const AxisAlignedBox3 &box,
-              float &position
-              );
+                                                        const AxisAlignedBox3 &box,
+                                                        float &position,
+                                                        int &raysBack,
+                                                        int &raysFront,
+                                                        int &pvsBack,
+                                                        int &pvsFront
+                                                        );
 
   void
@@ -590 +619 @@
   }
   
+        float
+        BestCostRatio(
+                                                                VssTreeLeaf *node,
+                                                                int &axis,
+                                                                float &position,
+                                                                int &raysBack,
+                                                                int &raysFront
+                                                                );
+        
+        float
+        EvalCostRatio(
+                                                                VssTreeLeaf *node,
+                                                                const int axis,
+                                                                const float position,
+                                                                int &raysBack,
+                                                                int &raysFront,
+                                                                int &pvsBack,
+                                                                int &pvsFront
+                                                                );
 
   AxisAlignedBox3 GetBBox(const VssTreeNode *node) {
@@ -664 +712 @@
   EvaluateLeafStats(const TraversalData &data);
 
+        void
+        EvalLeafPvs(VssTreeLeaf *leaf);
+
+        
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/default.env

@@ -9 +9 @@
 #       filename vienna.x3d
 # filename ../data/vienna/vienna-simple.x3d
-filename ../data/vienna/vienna-buildings.x3d
+# filename ../data/vienna/vienna-buildings.x3d
 # filename ../data/vienna/viewcells-25-sel.x3d
 # filename ../data/atlanta/atlanta2.x3d
 #       filename ../data/soda/soda.dat
-# filename ../data/soda/soda5.dat
+filename ../data/soda/soda5.dat
 }
 
@@ -22 +22 @@
 
 VssPreprocessor {
-        totalSamples 1000000
-        samplesPerPass  100000
+        totalSamples 200000
+        samplesPerPass  50000
+}
+
+VssTree {
+        epsilon         1e-6
+
+        maxDepth        40
+        minPvs          10
+        minSize         0.00001
+        maxCostRatio    1.0
+        maxRayContribution 0.5
+        
+        maxTotalMemory  500
+        maxStaticMemory 200
+
+        splitType regular
+#       splitType heuristics
+
+        numberOfEndPointDomains 10000
+        ct_div_ci       0.0
+        randomize       false
+
+        refDirBoxMaxSize        0.1
 }
 
@@ -66 +88 @@
 
 Sampling {
-        totalSamples 100000000
+        totalSamples 10000000
         samplesPerPass  3
 }

trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro

@@ -40 +40 @@
 ViewCell.cpp ViewCellBsp.cpp Halton.cpp VssRay.cpp VssTree.cpp VssPreprocessor.cpp
 
-