Changeset 448 for trunk/VUT/GtpVisibilityPreprocessor

Timestamp:

12/05/05 04:42:54 (19 years ago)

Author:

mattausch

Message:

fixed bug in VspBspTree?
view cells in VssPreprocessor?
bounding rays for vspkdtree

Location:

trunk/VUT/GtpVisibilityPreprocessor

Files:

• scripts/Preprocessor.vcproj (modified) (1 diff)
• scripts/default.env (modified) (8 diffs)
• src/AxisAlignedBox3.cpp (modified) (2 diffs)
• src/AxisAlignedBox3.h (modified) (1 diff)
• src/Environment.cpp (modified) (2 diffs)
• src/Exporter.h (modified) (5 diffs)
• src/Plane3.h (modified) (1 diff)
• src/Polygon3.cpp (modified) (13 diffs)
• src/Polygon3.h (modified) (4 diffs)
• src/Preprocessor.cpp (modified) (2 diffs)
• src/Ray.cpp (modified) (1 diff)
• src/Ray.h (modified) (2 diffs)
• src/RayInfo.cpp (modified) (2 diffs)
• src/RenderSimulator.cpp (modified) (4 diffs)
• src/RenderSimulator.h (modified) (2 diffs)
• src/RssPreprocessor.cpp (modified) (6 diffs)
• src/SamplingPreprocessor.cpp (modified) (5 diffs)
• src/Vector3.cpp (modified) (1 diff)
• src/Vector3.h (modified) (2 diffs)
• src/ViewCellBsp.cpp (modified) (51 diffs)
• src/ViewCellBsp.h (modified) (9 diffs)
• src/ViewCellsManager.cpp (modified) (34 diffs)
• src/ViewCellsManager.h (modified) (3 diffs)
• src/VspBspTree.cpp (modified) (5 diffs)
• src/VspBspTree.h (modified) (20 diffs)
• src/VspKdTree.cpp (modified) (4 diffs)
• src/VssPreprocessor.cpp (modified) (6 diffs)
• src/VssRay.h (modified) (1 diff)
• src/X3dExporter.h (modified) (5 diffs)
• src/main.cpp (modified) (1 diff)

trunk/VUT/GtpVisibilityPreprocessor/scripts/Preprocessor.vcproj

@@ -255 +255 @@
                         </File>
                         <File
+                                RelativePath="..\src\RssPreprocessor.cpp">
+                        </File>
+                        <File
+                                RelativePath="..\src\RssPreprocessor.h">
+                        </File>
+                        <File
+                                RelativePath="..\src\RssTree.cpp">
+                        </File>
+                        <File
+                                RelativePath="..\src\RssTree.h">
+                        </File>
+                        <File
                                 RelativePath="..\src\SamplingPreprocessor.cpp">
                         </File>

trunk/VUT/GtpVisibilityPreprocessor/scripts/default.env

@@ -19 +19 @@
 
 Preprocessor {
-        type sampling
-#       type vss
+#       type sampling
+        type vss
 }
 
@@ -94 +94 @@
 
 Sampling {
-        totalSamples 600000
+        totalSamples 300000
         samplesPerPass  3
 }
@@ -101 +101 @@
         loadFromFile false
         #type kdTree
-        #type vspTree
+        type vspKdTree
         #type bspTree
-        type vspBspTree
+        #type vspBspTree
         
         #type sceneDependent
@@ -131 +131 @@
 BspTree {
         Construction {
-                samples 200000
-                sideTolerance 0.005
+                samples 50000
+                epsilon 0.005
         }
 
@@ -193 +193 @@
                 minRays 200
                 minPolygons -1
-                maxDepth 50
+                maxDepth 40
                 minPvs 100
                 minArea 0.01
@@ -256 +256 @@
 
         Construction {
-                samples 500000
+                samples 100000
         }
         
@@ -278 +278 @@
 VspBspTree {
         Construction {
-                samples 200000
-                sideTolerance 0.005
+                samples 100000
+                epsilon 0.005
         }
 
@@ -309 +309 @@
                 minRays 200
                 minPolygons -1
-                maxDepth 50
+                maxDepth 40
                 minPvs 100
                 minArea 0.01

trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp

@@ -1662 +1662 @@
 
 // TODO: use a table to avoid normal and distance computations
-Polygon3 *AxisAlignedBox3::CrossSection(const Plane3 &plane)
+Polygon3 *AxisAlignedBox3::CrossSection(const Plane3 &plane) const
 {
         Polygon3 *planePoly = new Polygon3();
@@ -1746 +1746 @@
         return planePoly;
 }
+
+bool AxisAlignedBox3::GetRaySegment(const Ray &ray, 
+                                                                        float &minT, 
+                                                                        float &maxT) const
+{
+        maxT = 1e6;
+        minT = 0;
+        
+        // test with  bounding box
+        if (!GetMinMaxT(ray, &minT, &maxT))
+                return false;
+
+        if (minT < 0) // start ray from origin
+                minT = 0;
+
+        // bound ray or line segment
+        if (//(ray.GetType() == Ray::LOCAL_RAY) && 
+            !ray.intersections.empty() && 
+                (ray.intersections[0].mT <= maxT))
+        {
+                maxT = ray.intersections[0].mT;
+        }
+
+        return true;
+}

trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h

@@ -312 +312 @@
           @returns the cross section
   */
-  Polygon3 *CrossSection(const Plane3 &plane);
+  Polygon3 *CrossSection(const Plane3 &plane) const;
+
+  /** Computes minimal and maximal t of ray, including the object intersections.
+          @returns true if ray hits the bounding box.
+  */
+  bool GetRaySegment(const Ray &ray, 
+                                float &minT, 
+                                float &maxT) const;
 
 #define __EXTENT_HACK

trunk/VUT/GtpVisibilityPreprocessor/src/Environment.cpp

@@ -1269 +1269 @@
           "100000");
 
-  RegisterOption("BspTree.Construction.sideTolerance",
+  RegisterOption("BspTree.Construction.epsilon",
           optFloat,
           "-bsp_construction_side_tolerance=",
@@ -1732 +1732 @@
           "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");
-   RegisterOption("VspBspTree.Factor.pvs", optFloat, "-vsp_bsp_factor_pvs=", "1.0");
+  RegisterOption("VspBspTree.Construction.samples",
+          optInt,
+          "-bsp_construction_samples=",
+          "100000");
+
+  RegisterOption("VspBspTree.Construction.epsilon",
+          optFloat,
+          "-vsp_bsp_construction_side_tolerance=",
+          "0.002");
+
+  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");
+  RegisterOption("VspBspTree.Factor.pvs", optFloat, "-vsp_bsp_factor_pvs=", "1.0");
 
    //////////////////////////////////////////////////////////////////////////////////

trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h

@@ -20 +20 @@
 class Polygon3;
 class VssTree;
+class VspBspTree;
 class RssTree;
 
@@ -44 +45 @@
   virtual bool ExportScene(SceneGraphNode *node) = 0;
 
-
   virtual bool
   ExportBox(const AxisAlignedBox3 &box) = 0;
@@ -56 +56 @@
   virtual bool
   ExportVssTree2(const VssTree &tree,
-                                 const Vector3 direction
-                                 ) = 0;
-
-  virtual bool
-  ExportRssTree2(const RssTree &tree,
                                  const Vector3 direction
                                  ) = 0;
@@ -98 +93 @@
   virtual void 
   ExportBspSplits(const BspTree &tree, const bool exportDepth = false) = 0;
+
   virtual void
   ExportLeavesGeometry(const BspTree &tree, const vector<BspLeaf *> &leaves) = 0;
-        
+  
+  virtual void 
+  ExportBspSplits(const VspBspTree &tree, const bool exportDepth = false) = 0;
+
+  virtual void
+  ExportBspViewCellPartition(const VspBspTree &tree, const int maxPvs = 0) = 0;
+
   virtual void 
   ExportPolygons(const PolygonContainer &polys) = 0;
@@ -112 +114 @@
   virtual void
   ExportGeometry(const ObjectContainer &objects) = 0;
+
+  virtual bool
+  ExportRssTree2(const RssTree &tree,
+                                 const Vector3 direction) = 0;
 
   void SetExportRayDensity(const bool d) { mExportRayDensity = d; }

trunk/VUT/GtpVisibilityPreprocessor/src/Plane3.h

@@ -51 +51 @@
   */
   Vector3 FindIntersection(const Vector3 &a,
-                           const Vector3 &b,
-                           float *t = NULL,
-                           bool *coplanar = NULL
-                           ) const 
+                                                   const Vector3 &b,
+                                                   float *t = NULL,
+                                                   bool *coplanar = NULL) const 
   {
     const Vector3 v = b - a; // line from A to B

trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.cpp

@@ -42 +42 @@
 void Polygon3::Split(const Plane3 &partition, 
                                          Polygon3 &front, 
-                                         Polygon3 &back)
+                                         Polygon3 &back,
+                                         const float epsilon)
 {
         Vector3 ptA = mVertices.back();
         
-        int sideA = partition.Side(ptA, Vector3::sDistTolerance);
+        int sideA = partition.Side(ptA, epsilon);
         
         VertexContainer::const_iterator it;
@@ -53 +54 @@
 
         bool foundSplit = false;
-        // find line - plane intersections
+
+        //-- find line - plane intersections
         for (it = mVertices.begin(); it != mVertices.end(); ++ it)
         {
                 Vector3 ptB = *it;
-                int sideB = partition.Side(ptB, Vector3::sDistTolerance);
+                int sideB = partition.Side(ptB, epsilon);
         
                 // vertices on different sides => split
@@ -68 +70 @@
                         
                                 // test if split point not too close to previous split point
-                                if (!foundSplit || 
-                                        (SqrDistance(splitPt, lastSplit) > Vector3::sDistToleranceSqrt))
+                                if (!foundSplit || (!EpsilonEqualV3(splitPt, lastSplit, epsilon)))
                                 {
                                         // add vertex to both polygons
@@ -89 +90 @@
                                 // test if split point not too close to other split point
                                         // test if split point not too close to previous split point
-                                if (!foundSplit || 
-                                        (SqrDistance(splitPt, lastSplit) > Vector3::sDistToleranceSqrt))
+                                if (!foundSplit || (!EpsilonEqualV3(splitPt, lastSplit, epsilon)))
                                 {
                                         // add vertex to both polygons
@@ -125 +125 @@
 }
 
-int Polygon3::Side(const Plane3 &plane) const
-{
-        int classification = ClassifyPlane(plane);
+int Polygon3::Side(const Plane3 &plane, 
+                                   const float epsilon) const
+{
+        int classification = ClassifyPlane(plane, epsilon);
         
         if (classification == BACK_SIDE)
@@ -137 +138 @@
 }
 
-int Polygon3::ClassifyPlane(const Plane3 &plane, const float tolerance) const
+int Polygon3::ClassifyPlane(const Plane3 &plane, 
+                                                        const float epsilon) const
 {
         VertexContainer::const_iterator it;
@@ -148 +150 @@
         for (it = mVertices.begin(); it != mVertices.end(); ++ it)
         {
-                const int side = plane.Side(*it, tolerance);
+                const int side = plane.Side(*it, epsilon);
                 
                 if (side > 0)
@@ -179 +181 @@
 }
 
-int Polygon3::ClassifyPlane(const Plane3 &plane) const
-{
-        return ClassifyPlane(plane, Vector3::sDistTolerance);
-}
+/*int Polygon3::ClassifyPlane(const Plane3 &plane) const
+{
+        return ClassifyPlane(plane, Limits::Small);
+}*/
 
 
@@ -207 +209 @@
 }
 
-bool Polygon3::Valid() const
+bool Polygon3::Valid(const float epsilon) const
 {
         if (mVertices.size() < 3)
@@ -225 +227 @@
         for (it = mVertices.begin(); it != it_end; ++it)
         {
-                if (!(SqrDistance(vtx, *it) > sDistToleranceSqrt))
-                {
-                        Debug << "Malformed vertices:\n" << *this << endl;
+                if (!(EpsilonEqual(vtx, *it)))
+                {
+                        //Debug << "Malformed vertices:\n" << *this << endl;
                         return false;
                 }
@@ -385 +387 @@
 }
 
-int Polygon3::ClassifyPlane(const PolygonContainer &polys, const Plane3 &plane)
+int Polygon3::ClassifyPlane(const PolygonContainer &polys, 
+                                                        const Plane3 &plane,
+                                                        const float epsilon)
 {
         PolygonContainer::const_iterator it;
@@ -395 +399 @@
         for (it = polys.begin(); it != polys.end(); ++ it)
         {
-        const int cf = (*it)->ClassifyPlane(plane);
+        const int cf = (*it)->ClassifyPlane(plane, epsilon);
                 
         if (cf == FRONT_SIDE)
@@ -445 +449 @@
 
         for (pit = cell.begin(); pit != cell.end(); ++ pit)
+        {
+                if ((*pit)->mVertices.size() < 3)
+                        Debug << "ERROR" << endl;
+
                 area += (*pit)->GetArea();
+        }
 
         return area;

trunk/VUT/GtpVisibilityPreprocessor/src/Polygon3.h

@@ -50 +50 @@
                 @param front returns the front the front polygon
                 @param back returns the back polygon
+                @param epsilon epsilon where two points are considered equal
         */
         void Split(const Plane3 &partition, 
                            Polygon3 &front, 
-                           Polygon3 &back);
+                           Polygon3 &back,
+                           const float epsilon);// = Limits::Small);
 
         /** Returns the area of this polygon.
@@ -62 +64 @@
             @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
         */
-        int ClassifyPlane(const Plane3 &plane) const;
-
+        //int ClassifyPlane(const Plane3 &plane) const;
         
         /** Classify polygon with respect to the plane.
-            @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
+            @param epsilon tolerance value
+                @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
+
         */
-        int ClassifyPlane(const Plane3 &plane, const float tolerance) const;
+        int ClassifyPlane(const Plane3 &plane, const float epsilon) const;
 
         /** Side of the polygon with respect to the plane.
                 @returns 1 if on front side, -1 if on back side, 0 else.
         */
-        int Side(const Plane3 &plane) const;
+        int Side(const Plane3 &plane, const float epsilon) const;
 
         /** Scales the polygon about its center
@@ -85 +88 @@
                 @returns true if polygon is valid.
         */
-        bool Valid() const; 
+        bool Valid(const float epsilon) const; 
 
         /** Returns the surface normal.
@@ -127 +130 @@
             @returns one of BACK_SIDE, FRONT_SIDE, SPLIT, COINCIDENT
         */
-        static int ClassifyPlane(const PolygonContainer &polys, const Plane3 &plane);
+        static int ClassifyPlane(const PolygonContainer &polys, 
+                                                         const Plane3 &plane,
+                                                         const float epsilon);
 
 

trunk/VUT/GtpVisibilityPreprocessor/src/Preprocessor.cpp

@@ -178 +178 @@
                 mViewCellsManager = new KdViewCellsManager(mKdTree);
         }
-        if (strcmp(viewCellsStr, "bspTree") == 0)
+        else if (strcmp(viewCellsStr, "bspTree") == 0)
         {
                 mBspTree = new BspTree();
+
+                Debug << "view cell type: Bsp" << endl;
 
                 environment->GetIntValue("BspTree.Construction.samples", constructionSamples);
                 mViewCellsManager = new BspViewCellsManager(mBspTree, constructionSamples);
         }
-        if (strcmp(viewCellsStr, "vspBspTree") == 0)
+        else if (strcmp(viewCellsStr, "vspBspTree") == 0)
         {
                 mVspBspTree = new VspBspTree();
+
+                Debug << "view cell type: VspBsp" << endl;
 
                 environment->GetIntValue("VspBspTree.Construction.samples", constructionSamples);
@@ -204 +208 @@
                 mBspTree = new BspTree();
 
+                Debug << "view cell type: Bsp" << endl;
                 environment->GetIntValue("BspTree.Construction.samples", constructionSamples);
                 mViewCellsManager = new BspViewCellsManager(mBspTree, constructionSamples);

trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp

@@ -225 +225 @@
 mType(LOCAL_RAY)
 {
-        const float dist = Distance(vssRay.mTermination, vssRay.mOrigin);
-
-        if (dist)
-                dir = vssRay.GetDir() / dist;
+        float len = vssRay.Length();
+
+        if (!len)
+                len = Limits::Small;
+
+        dir = vssRay.GetDir() / len;
 
         if (vssRay.mTerminationObject)
-                intersections.push_back(Intersection(dist, vssRay.mTerminationObject, 0));
+                intersections.push_back(Intersection(len, vssRay.mTerminationObject, 0));
 
         Precompute();

trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h

@@ -15 +15 @@
 class VssRay;
 
+
 // -------------------------------------------------------------------
 // CRay class.  A ray is defined by a location and a direction.
@@ -48 +49 @@
                 
     bool operator<(
-                                                                         const Intersection &b) const {
+                                   const Intersection &b) const {
+
       return 
                                 mT

trunk/VUT/GtpVisibilityPreprocessor/src/RayInfo.cpp

@@ -45 +45 @@
 float RayInfo::ExtrapOrigin(const int axis) const 
 {
-        return mRay->GetOrigin(axis) + GetMinT()*mRay->GetDir(axis);
+        return mRay->GetOrigin(axis) + GetMinT() * mRay->GetDir(axis);
 }
                 
 float RayInfo::ExtrapTermination(const int axis) const 
 {
-        return mRay->GetOrigin(axis) + GetMaxT()*mRay->GetDir(axis);
+        return mRay->GetOrigin(axis) + GetMaxT() * mRay->GetDir(axis);
 }
 
 Vector3 RayInfo::ExtrapOrigin() const 
 {
-        return mRay->GetOrigin() * mRay->GetDir();
+        return mRay->GetOrigin() + GetMinT() * mRay->GetDir();
 }
                 
@@ -134 +134 @@
         t = splitPlane.FindT(mRay->GetOrigin(), mRay->GetTermination());
 
+//      Debug << "t: " << t << " mint " << GetMinT() << " maxt " << GetMaxT() << " orig: " << mRay->GetOrigin() <<
+//              " term " << mRay->GetTermination() << endl;
         if ((t < GetMinT()) || (t > GetMaxT()))
                 return splitPlane.Side(ExtrapOrigin());

trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.cpp

@@ -3 +3 @@
 #include "ViewCellBsp.h"
 #include "ViewCell.h"
+#include "VspBspTree.h"
+
 
 void SimulationStatistics::Print(ostream &app) const
@@ -53 +55 @@
 /*     class BspRenderSimulator implementation          */
 /********************************************************/
+
+
 BspRenderSimulator::BspRenderSimulator(BspTree *bspTree):
 mBspTree(bspTree)
@@ -139 +143 @@
 
 Real BspRenderSimulator::RenderPvs(ViewCell &viewCell, 
-                                                                           float objRenderTime) const
+                                                                   float objRenderTime) const
 {
         return viewCell.GetPvs().GetSize() * objRenderTime;
@@ -234 +238 @@
         return leaf->mKdPvs.GetSize() * objRenderTime;
 }
+
+/********************************************************/
+/*     class BspRenderSimulator implementation          */
+/********************************************************/
+VspBspRenderSimulator::VspBspRenderSimulator(VspBspTree *vspBspTree):
+mVspBspTree(vspBspTree)
+{
+}
+
+VspBspRenderSimulator::VspBspRenderSimulator(float objRenderCost, 
+                                                                                         float vcOverhead, 
+                                                                                         float moveSpeed,
+                                                                                         VspBspTree *vspBspTree):
+RenderSimulator(objRenderCost, vcOverhead, moveSpeed), 
+mVspBspTree(vspBspTree)
+{
+}
+
+SimulationStatistics VspBspRenderSimulator::SimulateRendering()
+{
+        SimulationStatistics simStat;
+
+        simStat.Reset();
+        simStat.Start();
+
+        Real renderTime = 0;
+        
+        // overhead for loading the PVS of the view cells
+        float loadPvsOverhead = 0;
+        // probability that view point lies in a view cell
+        float pInVcTotal = 0;
+
+        // total probability that a view cell border is crossed
+        const float pCrossVcTotal = mVspBspTree->GetBoundingBox().SurfaceArea();
+
+        // collect all view cells
+        ViewCellContainer viewCells;
+        mVspBspTree->CollectViewCells(viewCells);
+
+        ViewCellContainer::const_iterator it, it_end = viewCells.end();
+
+        // surface area substitute for probability
+        PolygonContainer geom;
+float overallarea = 0;
+        for (it = viewCells.begin(); it != it_end; ++ it)
+        {
+                // compute view cell area
+                mVspBspTree->ConstructGeometry(dynamic_cast<BspViewCell *>(*it), geom);
+
+                const float area = Polygon3::GetArea(geom);
+                if (area < 0.0001)
+                        Debug << "warning, area: " << area << endl;
+                CLEAR_CONTAINER(geom);
+
+                // area substitute for view point probability
+                float pInVc = area;
+                                
+                // compute render time of PVS times probability that view point is in view cell
+                float vcCost = pInVc * RenderPvs(*(*it), mObjRenderCost);
+                //Debug << "p: " << pInVc << " rendercost: " << RenderPvs(*(*it), mObjRenderCost) << endl;
+                renderTime += vcCost;
+
+                if (vcCost > simStat.maxCost)
+                        simStat.maxCost = vcCost;
+                else if (vcCost < simStat.minCost)
+                        simStat.minCost = vcCost;
+
+                // probability that a view cell border is crossed
+                float pCrossVc = area * mMoveSpeed;
+
+                // crossing the border of a view cell is also depending on the move speed
+                loadPvsOverhead += pCrossVc * mVcOverhead;
+overallarea += area;
+                pInVcTotal += pInVc;
+        }
+        Debug << "overall area: " << overallarea << endl;
+        
+        renderTime /= pInVcTotal;
+        loadPvsOverhead /= pCrossVcTotal;
+
+        simStat.avgRtWithoutOverhead = renderTime;
+        simStat.avgRenderTime = renderTime + loadPvsOverhead;
+        
+        simStat.maxCost /= pInVcTotal;
+        simStat.minCost /= pInVcTotal;
+
+        simStat.Stop();
+
+        return simStat;
+}
+
+Real VspBspRenderSimulator::RenderPvs(ViewCell &viewCell, 
+                                                                          float objRenderTime) const
+{
+        return viewCell.GetPvs().GetSize() * objRenderTime;
+}

trunk/VUT/GtpVisibilityPreprocessor/src/RenderSimulator.h

@@ -9 +9 @@
 class ViewCell;
 class KdLeaf;
+class VspBspTree;
 
 class SimulationStatistics: public StatisticsBase
@@ -105 +106 @@
         Real RenderPvs(ViewCell &viewCell, const float objRenderTime) const;
 
+private:
         BspTree *mBspTree;
+};
+
+class VspBspRenderSimulator: public RenderSimulator
+{
+public:
+        VspBspRenderSimulator(VspBspTree *vspBspTree);
+
+        VspBspRenderSimulator(float objRenderCost, 
+                                                  float vcOverhead, 
+                                                  float moveSpeed, 
+                                                  VspBspTree *vspBspTree);
+
+        SimulationStatistics SimulateRendering();
+
+protected:
+        /** Simulates rendering of the pvs of one view cell, with given rendering time for an object.
+                @param viewCell the view cell holding the Pvs
+                @param objRenderTime estimated render time for one object of the Pvs
+        */
+        Real RenderPvs(ViewCell &viewCell, const float objRenderTime) const;
+
+private:
+        VspBspTree *mVspBspTree;
 };
 

trunk/VUT/GtpVisibilityPreprocessor/src/RssPreprocessor.cpp

@@ -9 +9 @@
 #include "VssRay.h"
 #include "RssTree.h"
-
+#include "ViewCellBsp.h"
 
 static bool useViewSpaceBox = true;//true;
@@ -474 +474 @@
   }
   
-  //-- construct BSP view cells
-  if (ViewCell::sHierarchy == ViewCell::BSP) 
-  {
-        const int bspSamples = min((int)mVssRays.size(), mBspConstructionSamples);
-        
-        Debug << "bpssamples: " << bspSamples << endl;
-        for (int i = 0; i < bspSamples; ++ i)
-          bspRays.push_back(new Ray(*mVssRays[i]));
-        
-        //-- construct BSP tree using the samples
-        mBspTree = new BspTree(&mUnbounded);    
-        
-        ObjectContainer objects;
-        mSceneGraph->CollectObjects(&objects);
-        mBspTree->Construct(objects, bspRays);
-  }
-  
   rssTree = new RssTree;
   // viewcells = Construct(mVssRays);
@@ -567 +550 @@
         }
 
-        // cast rays into BSP tree
-        if (ViewCell::sHierarchy == ViewCell::BSP) 
-          {
-                for (int i = 0; i < (int)vssRays.size(); ++ i)
-                  {
-                        CastRay(*mBspTree, *vssRays[i]);
-                  }
-          }
-
         samples+=num;
         float pvs = rssTree->GetAvgPvsSize();
@@ -598 +572 @@
   delete rssTree;
   
-  if (ViewCell::sHierarchy == ViewCell::BSP) 
-        {
-          Debug << mBspTree->GetStatistics();
-
-          ObjectContainer objects;
-          ExportSplits(objects, bspRays, 10000);
-          ExportBspPvs(objects, bspRays, 10000);
-
-          BspViewCellsStatistics stat;
-          mBspTree->EvaluateViewCellsStats(stat);
-          Debug << "original view cell partition:\n" << stat << endl;
-
-          // clear BSP samples
-          CLEAR_CONTAINER(bspRays);
-        }
 
   return true;
@@ -646 +605 @@
           BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
           // if ray not in unbounded space
-          if (leaf->GetViewCell() != &mUnbounded)
+/*        if (leaf->GetViewCell() != &mUnbounded)
                 contributingSamples += 
-                  leaf->GetViewCell()->GetPvs().AddSample(obj);
+                  leaf->GetViewCell()->GetPvs().AddSample(obj);*/
         }
  
@@ -657 +616 @@
                 BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
                 
-                if (leaf->GetViewCell() != &mUnbounded)
+        /*      if (leaf->GetViewCell() != &mUnbounded)
                   leaf->GetViewCell()->
-                        AddPassingRay(ray, contributingSamples ? 1 : 0);
+                        AddPassingRay(ray, contributingSamples ? 1 : 0);*/
           }
   

trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp

@@ -75 +75 @@
 
         mViewCellsManager->ComputeSampleContributions(rays, 
-                                                                                                  castRaysToViewCells,
                                                                                                   sampleContributions, 
-                                                                                                  contributingSamples);
+                                                                                                  contributingSamples,
+                                                                                                  castRaysToViewCells);
 }
 
@@ -397 +397 @@
         //-- several visualizations and statistics
                  
-        //-- render simulation
-        cout << "\nevaluating bsp view cells render time before merge ... ";
+        //-- render simulation after merge
+        cout << "\nevaluating bsp view cells render time after merge ... ";
                  
         const SimulationStatistics ss = mViewCellsManager->SimulateRendering();
@@ -405 +405 @@
         cout << ss << endl;
         Debug << ss << endl;
-                 
+        
+
         mViewCellsManager->Visualize(objects, mSampleRays);     
     
@@ -424 +425 @@
                 RayContainer::const_iterator it, it_end = newRays.end();
 
-                if ((int)mVssSampleRays.size() < mViewCellsManager->GetConstructionSamples())
+                if ((int)mVssSampleRays.size() < 
+                        mViewCellsManager->GetConstructionSamples())
                 {
                         for (it = newRays.begin(); it != it_end; ++ it)
@@ -433 +435 @@
                 else
                 {
-                        // construct view cells using the collected samples
-                        cout << "building view cells from " << (int)mVssSampleRays.size() << " samples " << endl;
-
+                        // construct view cells
                         mViewCellsManager->Construct(objects, mVssSampleRays);
-                        
+                
                         // throw away samples 
-                        CLEAR_CONTAINER(mVssSampleRays);
-                }
-        }
-        // Need (ordered) rays for post processing => collect new rays
+                        //CLEAR_CONTAINER(mVssSampleRays);
+                }
+        }
+        // Need rays (with ordered intersections) for post processing => collect new rays
         else if (((int)mSampleRays.size() < mViewCellsManager->GetPostProcessSamples()) ||
                          ((int)mSampleRays.size() < mViewCellsManager->GetVisualizationSamples()))

trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.cpp

@@ -2 +2 @@
 #include "Vector3.h"
 #include "Halton.h"
-
-float Vector3::sDistTolerance = 0.002f;
-float Vector3::sDistToleranceSqrt = 0.000004f;
 
 // Given min a vector to minimize and a candidate vector, replace

trunk/VUT/GtpVisibilityPreprocessor/src/Vector3.h

@@ -36 +36 @@
   Vector3(float X) { x = y = z = X; }
   Vector3(const Vector3 &v) { x = v.x; y = v.y; z = v.z; }
-
-  /// the distance where two points are still considered equal
-  static float sDistTolerance;
-  static float sDistToleranceSqrt;
 
   // Functions to get at the vector components
@@ -480 +476 @@
 EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2)
 {
-  return EpsilonEqualV3(v1,v2,Limits::Small);
+  return EpsilonEqualV3(v1, v2, Limits::Small);
 }
 

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp

@@ -98 +98 @@
 }
 
-Plane3 *BspInterior::GetPlane()
-{
-        return &mPlane;
+Plane3 BspInterior::GetPlane() const
+{
+        return mPlane;
 }
 
@@ -117 +117 @@
 }
 
-int BspInterior::SplitPolygons(PolygonContainer &polys, 
-                                                           PolygonContainer &frontPolys, 
-                                                           PolygonContainer &backPolys, 
-                                                           PolygonContainer &coincident)
-{
-        int splits = 0;
-
-#ifdef _Debug
-        Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
-#endif
-        while (!polys.empty())
-        {
-                Polygon3 *poly = polys.back();
-                polys.pop_back();
-
-                //Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";
-
-                // classify polygon
-                const int cf = poly->ClassifyPlane(mPlane);
-
-                switch (cf)
-                {
-                        case Polygon3::COINCIDENT:
-                                coincident.push_back(poly);
-                                break;                  
-                        case Polygon3::FRONT_SIDE:      
-                                frontPolys.push_back(poly);
-                                break;
-                        case Polygon3::BACK_SIDE:
-                                backPolys.push_back(poly);
-                                break;
-                        case Polygon3::SPLIT:
-                                {
-                                        Polygon3 *front_piece = new Polygon3(poly->mParent);
-                                        Polygon3 *back_piece = new Polygon3(poly->mParent);
-
-                                        //-- split polygon into front and back part
-                                        poly->Split(mPlane, *front_piece, *back_piece);
-                                        
-                                        ++ splits; // increase number of splits
-
-                                        //-- inherit rays from parent polygon for blocked ray criterium
-                                        poly->InheritRays(*front_piece, *back_piece);
-                                        //Debug << "p: " << poly->mPiercingRays.size() << " f: " << front_piece->mPiercingRays.size() << " b: " << back_piece->mPiercingRays.size() << endl;
-                                
-                                        // check if polygons still valid 
-                                        if (front_piece->Valid())
-                                                frontPolys.push_back(front_piece);
-                                        else
-                                                DEL_PTR(front_piece);
-                                
-                                        if (back_piece->Valid())
-                                                backPolys.push_back(back_piece);
-                                        else                            
-                                                DEL_PTR(back_piece);
-                                
-#ifdef _DEBUG
-                                        Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
-#endif
-                                        DEL_PTR(poly);
-                                }
-                                break;
-                        default:
-                Debug << "SHOULD NEVER COME HERE\n";
-                                break;
-                }
-        }
-
-        return splits;
-}
-
 /****************************************************************/
 /*                  class BspLeaf implementation                */
 /****************************************************************/
+
+
 BspLeaf::BspLeaf(): mViewCell(NULL)
 {
 }
+
 
 BspLeaf::BspLeaf(BspViewCell *viewCell): 
@@ -199 +131 @@
 {
 }
+
 
 BspLeaf::BspLeaf(BspInterior *parent): 
@@ -204 +137 @@
 {}
 
+
+
 BspLeaf::BspLeaf(BspInterior *parent, BspViewCell *viewCell): 
 BspNode(parent), mViewCell(viewCell)
@@ -222 +157 @@
 { 
         return true; 
-}
-
-void BspLeaf::AddToPvs(const BoundedRayContainer &rays, 
-                                           int &sampleContributions,
-                                           int &contributingSamples)
-{
-        sampleContributions = 0;
-        contributingSamples = 0;
-
-    BoundedRayContainer::const_iterator it, it_end = rays.end();
-
-        // add contributions from samples to the PVS
-        for (it = rays.begin(); it != it_end; ++ it)
-        {
-                int contribution = 0;
-                Ray *ray = (*it)->mRay;
-                        
-                if (!ray->intersections.empty())
-                        contribution += mViewCell->GetPvs().AddSample(ray->intersections[0].mObject);
-                
-                if (ray->sourceObject.mObject)
-                        contribution += mViewCell->GetPvs().AddSample(ray->sourceObject.mObject);
-
-                if (contribution)
-                {
-                        sampleContributions += contribution;
-                        ++ contributingSamples;
-                }
-
-                //if (ray->mFlags & Ray::STORE_BSP_INTERSECTIONS)
-                        ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
-        }
 }
 
@@ -308 +211 @@
         environment->GetIntValue("BspTree.maxTests", mMaxTests);
 
-        environment->GetFloatValue("BspTree.Construction.sideTolerance", Vector3::sDistTolerance);
-        Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
-
-
+        environment->GetFloatValue("BspTree.Construction.epsilon", mEpsilon);
+        
     Debug << "BSP max depth: " << mTermMaxDepth << endl;
         Debug << "BSP min PVS: " << mTermMinPvs << endl;
@@ -347 +248 @@
 }
 
+
 const BspTreeStatistics &BspTree::GetStatistics() const 
 {
         return mStat;
 }
+
+
+int BspTree::SplitPolygons(const Plane3 &plane,
+                                                   PolygonContainer &polys, 
+                                                   PolygonContainer &frontPolys, 
+                                                   PolygonContainer &backPolys, 
+                                                   PolygonContainer &coincident) const
+{
+        int splits = 0;
+
+#ifdef _Debug
+        Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
+#endif
+        while (!polys.empty())
+        {
+                Polygon3 *poly = polys.back();
+                polys.pop_back();
+
+                //Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";
+
+                // classify polygon
+                const int cf = poly->ClassifyPlane(plane, mEpsilon);
+
+                switch (cf)
+                {
+                        case Polygon3::COINCIDENT:
+                                coincident.push_back(poly);
+                                break;                  
+                        case Polygon3::FRONT_SIDE:      
+                                frontPolys.push_back(poly);
+                                break;
+                        case Polygon3::BACK_SIDE:
+                                backPolys.push_back(poly);
+                                break;
+                        case Polygon3::SPLIT:
+                                {
+                                        Polygon3 *front_piece = new Polygon3(poly->mParent);
+                                        Polygon3 *back_piece = new Polygon3(poly->mParent);
+
+                                        //-- split polygon into front and back part
+                                        poly->Split(plane, 
+                                                                *front_piece, 
+                                                                *back_piece,
+                                                                mEpsilon);
+                                        
+                                        ++ splits; // increase number of splits
+
+                                        //-- inherit rays from parent polygon for blocked ray criterium
+                                        poly->InheritRays(*front_piece, *back_piece);
+                                
+                                        // check if polygons still valid 
+                                        if (front_piece->Valid(mEpsilon))
+                                                frontPolys.push_back(front_piece);
+                                        else
+                                                DEL_PTR(front_piece);
+                                
+                                        if (back_piece->Valid(mEpsilon))
+                                                backPolys.push_back(back_piece);
+                                        else                            
+                                                DEL_PTR(back_piece);
+                                
+#ifdef _DEBUG
+                                        Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
+#endif
+                                        DEL_PTR(poly);
+                                }
+                                break;
+                        default:
+                Debug << "SHOULD NEVER COME HERE\n";
+                                break;
+                }
+        }
+
+        return splits;
+}
+
 
 void BspTreeStatistics::Print(ostream &app) const
@@ -458 +436 @@
                 
                                 // split viewcell polygons with respect to split plane
-                                splits += interior->SplitPolygons(*tData.mPolygons,
-                                                                                                  *frontPolys, 
-                                                                                                  *backPolys,
-                                                                                                  coincident);
+                                splits += SplitPolygons(interior->GetPlane(),
+                                                                                *tData.mPolygons, 
+                                                                                *frontPolys, 
+                                                                                *backPolys, 
+                                                                                coincident);
                                 
                                 // extract view cells associated with the split polygons
@@ -505 +484 @@
                         // cleanup
                         DEL_PTR(tData.mPolygons);
+                        DEL_PTR(tData.mRays);
                 }
                 else
@@ -514 +494 @@
 }
 
-int BspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
+int BspTree::AddMeshToPolygons(Mesh *mesh, 
+                                                           PolygonContainer &polys, 
+                                                           MeshInstance *parent)
 {
         FaceContainer::const_iterator fi;
@@ -523 +505 @@
                 Polygon3 *poly = new Polygon3((*fi), mesh);
                 
-                if (poly->Valid())
+                if (poly->Valid(mEpsilon))
                 {
                         poly->mParent = parent; // set parent intersectable
@@ -684 +666 @@
 
                 float minT, maxT;
-                if (BoundRay(*ray, minT, maxT))
+                if (mBox.GetRaySegment(*ray, minT, maxT))
                         rays->push_back(new BoundedRay(ray, minT, maxT));
         }
@@ -719 +701 @@
 
                 float minT, maxT;
-                if (BoundRay(*ray, minT, maxT))
+                if (mBox.GetRaySegment(*ray, minT, maxT))
                         rays->push_back(new BoundedRay(ray, minT, maxT));
         }
@@ -801 +783 @@
                 {
                         int conSamp = 0, sampCon = 0;
-                        leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
+                        AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
                         
                         mStat.contributingSamples += conSamp;
@@ -818 +800 @@
                 DEL_PTR(tData.mPolygons);
                 DEL_PTR(tData.mRays);
-                        
+                DEL_PTR(tData.mGeometry);
+
                 return leaf;
         }
@@ -860 +843 @@
         // cleanup
         DEL_PTR(tData.mNode);
+
         DEL_PTR(tData.mPolygons);
         DEL_PTR(tData.mRays);
@@ -916 +900 @@
         
 
-        Debug << "number of rays: " << tData.mRays->size() << endl;
-        Debug << "number of polys: " << tData.mPolygons->size() << endl;
+        Debug << "number of rays: " << (int)tData.mRays->size() << endl;
+        Debug << "number of polys: " << (int)tData.mPolygons->size() << endl;
 
         startTime = GetTime();
@@ -928 +912 @@
         startTime = GetTime();
         // subdivide polygons with plane
-        mStat.splits += interior->SplitPolygons(*tData.mPolygons, 
-                                                    *frontData.mPolygons, 
-                                                                                        *backData.mPolygons, 
-                                                                                        coincident);
+        mStat.splits += SplitPolygons(interior->GetPlane(),
+                                                                  *tData.mPolygons, 
+                                          *frontData.mPolygons, 
+                                                                  *backData.mPolygons, 
+                                                                  coincident);
 
         Debug << "time used for polygon splitting: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
@@ -944 +929 @@
                 tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
                                                                            *backData.mGeometry, 
-                                                                           *this, 
-                                                                           interior->mPlane);
+                                                                           interior->mPlane,
+                                                                           mBox,
+                                                                           mEpsilon);
         
                 
@@ -974 +960 @@
         interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
         
-        frontData.mNode = interior->mFront;
-        backData.mNode = interior->mBack;
+        frontData.mNode = interior->GetFront();
+        backData.mNode = interior->GetBack();
         
         //DEL_PTR(leaf);
@@ -1269 +1255 @@
 
                         BoundedRay *bRay = (*rays)[idx[j]];
-                        pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : bRay->mRay->Extrap(bRay->mMaxT);
+                        pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : 
+                                                                bRay->mRay->Extrap(bRay->mMaxT);
                 }       
                         
@@ -1347 +1334 @@
                 Polygon3 *poly = polys[testIdx];
 
-        const int classification = poly->ClassifyPlane(candidatePlane);
+        const int classification = 
+                        poly->ClassifyPlane(candidatePlane, mEpsilon);
 
                 if (mSplitPlaneStrategy & BALANCED_POLYS)
@@ -1437 +1425 @@
 }
 
-bool BspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
-{
-        maxT = 1e6;
-        minT = 0;
-        
-        // test with tree bounding box
-        if (!mBox.GetMinMaxT(ray, &minT, &maxT))
-                return false;
-
-        if (minT < 0) // start ray from origin
-                minT = 0;
-
-        // bound ray or line segment
-        if (//(ray.GetType() == Ray::LOCAL_RAY) && 
-            !ray.intersections.empty() && 
-                (ray.intersections[0].mT <= maxT))
-        {
-                maxT = ray.intersections[0].mT;
-        }
-
-        return true;
-}
 
 inline void BspTree::GenerateUniqueIdsForPvs()
@@ -1499 +1465 @@
                         BspNodeGeometry backCell;
 
-                        cell.SplitGeometry(frontCell, backCell, *this, candidatePlane);
+                        cell.SplitGeometry(frontCell, 
+                                                           backCell, 
+                                                           candidatePlane,
+                                                           mBox,
+                                                           mEpsilon);
                 
                         pFront = frontCell.GetArea();
@@ -1803 +1773 @@
         float maxt, mint;
 
-        if (!BoundRay(ray, mint, maxt))
+        if (!mBox.GetRaySegment(ray, mint, maxt))
                 return 0;
 
@@ -1818 +1788 @@
                 if (!node->IsLeaf()) 
                 {
-                        BspInterior *in = (BspInterior *) node;
+                        BspInterior *in = dynamic_cast<BspInterior *>(node);
                         
-                        Plane3 *splitPlane = in->GetPlane();
-
-                        int entSide = splitPlane->Side(entp);
-                        int extSide = splitPlane->Side(extp);
-
-                        Vector3 intersection;
+                        Plane3 splitPlane = in->GetPlane();
+                        const int entSide = splitPlane.Side(entp);
+                        const int extSide = splitPlane.Side(extp);
 
                         if (entSide < 0)
@@ -1858 +1825 @@
                         // find intersection of ray segment with plane
                         float t;
-                        extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &t);
+                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
                         maxt *= t;
                         
@@ -1934 +1901 @@
                         BspInterior *interior = dynamic_cast<BspInterior *>(node);
 
-                        nodeStack.push(interior->mFront);
-                        nodeStack.push(interior->mBack);
+                        nodeStack.push(interior->GetFront());
+                        nodeStack.push(interior->GetBack());
                 }
         }
@@ -1989 +1956 @@
                         BspInterior *interior = dynamic_cast<BspInterior *>(node);
 
-                        nodeStack.push(interior->mFront);
-                        nodeStack.push(interior->mBack);
+                        nodeStack.push(interior->GetFront());
+                        nodeStack.push(interior->GetBack());
                 }
         }
@@ -2101 +2068 @@
                 {
                         BspInterior *interior = dynamic_cast<BspInterior *>(n);
-                        Plane3 halfSpace = *dynamic_cast<BspInterior *>(interior)->GetPlane();
-
-            if (interior->mFront != lastNode)
+                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
+
+            if (interior->GetFront() != lastNode)
                                 halfSpace.ReverseOrientation();
 
@@ -2135 +2102 @@
         ExtractHalfSpaces(n, halfSpaces);
 
-        PolygonContainer candidatePolys;
+        PolygonContainer candidates;
 
         // bounded planes are added to the polygons (reverse polygons
@@ -2143 +2110 @@
                 Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
                 
-                if (p->Valid())
-                {
-                        candidatePolys.push_back(p->CreateReversePolygon());
+                if (p->Valid(mEpsilon))
+                {
+                        candidates.push_back(p->CreateReversePolygon());
                         DEL_PTR(p);
                 }
@@ -2158 +2125 @@
                         vertices.push_back(mBox.GetFace(i).mVertices[j]);
 
-                candidatePolys.push_back(new Polygon3(vertices));
-        }
-
-        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
+                candidates.push_back(new Polygon3(vertices));
+        }
+
+        for (int i = 0; i < (int)candidates.size(); ++ i)
         {
                 // polygon is split by all other planes
-                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
+                for (int j = 0; (j < (int)halfSpaces.size()) && candidates[i]; ++ j)
                 {
                         if (i == j) // polygon and plane are coincident
@@ -2172 +2139 @@
                         Polygon3 *frontPoly, *backPoly;
 
-                        const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
+                        const int cf = candidates[i]->
+                                ClassifyPlane(halfSpaces[j], mEpsilon);
                         
                         switch (cf)
@@ -2180 +2148 @@
                                         backPoly = new Polygon3();
 
-                                        candidatePolys[i]->Split(halfSpaces[j], 
-                                                                                         *frontPoly, 
-                                                                                         *backPoly);
-
-                                        DEL_PTR(candidatePolys[i]);
-
-                                        if (frontPoly->Valid())
-                                                candidatePolys[i] = frontPoly;
+                                        candidates[i]->Split(halfSpaces[j], 
+                                                                                 *frontPoly, 
+                                                                                 *backPoly,
+                                                                                 mEpsilon);
+
+                                        DEL_PTR(candidates[i]);
+
+                                        if (frontPoly->Valid(mEpsilon))
+                                                candidates[i] = frontPoly;
                                         else
                                                 DEL_PTR(frontPoly);
@@ -2194 +2163 @@
                                         break;
                                 case Polygon3::BACK_SIDE:
-                                        DEL_PTR(candidatePolys[i]);
+                                        DEL_PTR(candidates[i]);
                                         break;
                                 // just take polygon as it is
@@ -2204 +2173 @@
                 }
                 
-                if (candidatePolys[i])
-                        cell.push_back(candidatePolys[i]);
+                if (candidates[i])
+                        cell.push_back(candidates[i]);
         }
 }
@@ -2242 +2211 @@
                                 {
                                         const int cf = 
-                                                Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);
+                                                Polygon3::ClassifyPlane(neighborCandidate, 
+                                                                                                halfSpaces[i],
+                                                                                                mEpsilon);
 
                                         if (cf == Polygon3::BACK_SIDE)
@@ -2258 +2229 @@
                         BspInterior *interior = dynamic_cast<BspInterior *>(node);
         
-                        const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);
+                        const int cf = Polygon3::ClassifyPlane(cell, 
+                                                                                                   interior->mPlane,
+                                                                                                   mEpsilon);
 
                         if (cf == Polygon3::FRONT_SIDE)
-                                nodeStack.push(interior->mFront);
+                                nodeStack.push(interior->GetFront());
                         else
                                 if (cf == Polygon3::BACK_SIDE)
-                                        nodeStack.push(interior->mBack);
+                                        nodeStack.push(interior->GetBack());
                                 else 
                                 {
                                         // random decision
-                                        nodeStack.push(interior->mBack);
-                                        nodeStack.push(interior->mFront);
+                                        nodeStack.push(interior->GetBack());
+                                        nodeStack.push(interior->GetFront());
                                 }
                 }
@@ -2305 +2278 @@
                         ConstructGeometry(interior, cell);
 
-                        const int cf = Polygon3::ClassifyPlane(cell, halfspace);
+                        const int cf = Polygon3::ClassifyPlane(cell, 
+                                                                                                   halfspace, 
+                                                                                                   mEpsilon);
 
                         if (cf == Polygon3::BACK_SIDE)
-                                next = interior->mFront;
+                                next = interior->GetFront();
                         else
                                 if (cf == Polygon3::FRONT_SIDE)
-                                        next = interior->mFront;
+                                        next = interior->GetFront();
                         else 
                         {
                                 // random decision
                                 if (mask & 1)
-                                        next = interior->mBack;
+                                        next = interior->GetBack();
                                 else
-                                        next = interior->mFront;
+                                        next = interior->GetFront();
                                 mask = mask >> 1;
                         }
@@ -2353 +2328 @@
                         // random decision
                         if (mask & 1)
-                                nodeStack.push(interior->mBack);
+                                nodeStack.push(interior->GetBack());
                         else
-                                nodeStack.push(interior->mFront);
+                                nodeStack.push(interior->GetFront());
 
                         mask = mask >> 1;
@@ -2362 +2337 @@
         
         return NULL;
+}
+
+void BspTree::AddToPvs(BspLeaf *leaf,
+                                           const BoundedRayContainer &rays, 
+                                           int &sampleContributions,
+                                           int &contributingSamples)
+{
+        sampleContributions = 0;
+        contributingSamples = 0;
+
+    BoundedRayContainer::const_iterator it, it_end = rays.end();
+
+        ViewCell *vc = leaf->GetViewCell();
+
+        // add contributions from samples to the PVS
+        for (it = rays.begin(); it != it_end; ++ it)
+        {
+                int contribution = 0;
+                Ray *ray = (*it)->mRay;
+                        
+                if (!ray->intersections.empty())
+                        contribution += vc->GetPvs().AddSample(ray->intersections[0].mObject);
+                
+                if (ray->sourceObject.mObject)
+                        contribution += vc->GetPvs().AddSample(ray->sourceObject.mObject);
+
+                if (contribution)
+                {
+                        sampleContributions += contribution;
+                        ++ contributingSamples;
+                }
+
+                //if (ray->mFlags & Ray::STORE_BSP_INTERSECTIONS)
+                //      ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
+        }
 }
 
@@ -2397 +2407 @@
 }
 
-/*************************************************************
- *            BspNodeGeometry Implementation                 *
- *************************************************************/
+float BspTree::GetEpsilon() const
+{
+        return mEpsilon;
+}
+
+void BspViewCellsStatistics::Print(ostream &app) const
+{
+        app << "===== BspViewCells statistics ===============\n";
+
+        app << setprecision(4);
+
+        //app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
+
+        app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
+
+        app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
+
+        app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
+
+        app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
+
+        app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
+
+        app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
+
+        app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgBspLeaves() << endl;
+
+        app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxBspLeaves << endl;
+        
+        app << "===== END OF BspViewCells statistics ==========\n";
+}
+
+/*************************************************************/
+/*            BspNodeGeometry Implementation                 */
+/*************************************************************/
 
 BspNodeGeometry::~BspNodeGeometry()
@@ -2413 +2455 @@
 void BspNodeGeometry::SplitGeometry(BspNodeGeometry &front,
                                                                         BspNodeGeometry &back,
-                                                                        const BspTree &tree,                                             
-                                                                        const Plane3 &splitPlane) const
+                                                                        const Plane3 &splitPlane,
+                                                                        const AxisAlignedBox3 &box,
+                                                                        const float epsilon) const
 {       
         // get cross section of new polygon
-        Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);
-
-        planePoly = SplitPolygon(planePoly, tree);
-
-        //-- plane poly splits all other cell polygons
+        Polygon3 *planePoly = box.CrossSection(splitPlane);
+
+        // split polygon with all other polygons
+        planePoly = SplitPolygon(planePoly, epsilon);
+
+        //-- new polygon splits all other polygons
         for (int i = 0; i < (int)mPolys.size(); ++ i)
         {
-                const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
+                /// don't use epsilon here to get exact split planes
+                const int cf = 
+                        mPolys[i]->ClassifyPlane(splitPlane, Limits::Small);
                         
-                // split new polygon with all previous planes
                 switch (cf)
                 {
@@ -2436 +2481 @@
                                         Polygon3 *backPoly = new Polygon3();
                                 
-                                        poly->Split(splitPlane, *frontPoly, *backPoly);
+                                        poly->Split(splitPlane, 
+                                                                *frontPoly, 
+                                                                *backPoly,
+                                                                epsilon);
 
                                         DEL_PTR(poly);
 
-                                        if (frontPoly->Valid())
+                                        if (frontPoly->Valid(epsilon))
                                                 front.mPolys.push_back(frontPoly);
                                         else
                                                 DEL_PTR(frontPoly);
 
-                                        if (backPoly->Valid())
+                                        if (backPoly->Valid(epsilon))
                                                 back.mPolys.push_back(backPoly);
                                         else
@@ -2467 +2515 @@
         }
 
-        //-- finally add the new polygon to the child cells
+        //-- finally add the new polygon to the child node geometries
         if (planePoly)
         {
@@ -2481 +2529 @@
 
 Polygon3 *BspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
-                                                                                const BspTree &tree) const
-{
+                                                                                const float epsilon) const
+{
+        if (!planePoly->Valid(epsilon))
+                DEL_PTR(planePoly);
+
         // polygon is split by all other planes
         for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
@@ -2488 +2539 @@
                 Plane3 plane = mPolys[i]->GetSupportingPlane();
 
+                /// don't use epsilon here to get exact split planes
                 const int cf = 
-                        planePoly->ClassifyPlane(plane, 0.00001f);
+                        planePoly->ClassifyPlane(plane, Limits::Small);
                         
                 // split new polygon with all previous planes
@@ -2499 +2551 @@
                                         Polygon3 *backPoly = new Polygon3();
 
-                                        planePoly->Split(plane, *frontPoly, *backPoly);
+                                        planePoly->Split(plane, 
+                                                                         *frontPoly, 
+                                                                         *backPoly,
+                                                                         epsilon);
                                         
                                         // don't need anymore
@@ -2506 +2561 @@
 
                                         // back polygon is belonging to geometry
-                                        if (backPoly->Valid())
+                                        if (backPoly->Valid(epsilon))
                                                 planePoly = backPoly;
                                         else
@@ -2525 +2580 @@
         return planePoly;
 }
-
-void BspViewCellsStatistics::Print(ostream &app) const
-{
-        app << "===== BspViewCells statistics ===============\n";
-
-        app << setprecision(4);
-
-        //app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
-
-        app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
-
-        app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
-
-        app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
-
-        app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
-
-        app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
-
-        app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
-
-        app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgBspLeaves() << endl;
-
-        app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxBspLeaves << endl;
-        
-        app << "===== END OF BspViewCells statistics ==========\n";
-}

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.h

@@ -7 +7 @@
 #include <stack>
 #include "Statistics.h"
+#include "VssRay.h"
+
 
 class ViewCell;
@@ -27 +29 @@
         float GetArea() const;
         
-        /** Computes new cell based on the old cell definition and a new split plane
-                @param side indicates which side of the halfspace 
-        */
-        void SplitGeometry(BspNodeGeometry &front, 
+        /** Computes new front and back geometry based on the old cell 
+                geometry and a new split plane
+        */
+        void SplitGeometry(BspNodeGeometry &front,
                                            BspNodeGeometry &back,
-                                           const BspTree &tree,
-                                           const Plane3 &splitPlane) const;
-
-        Polygon3 *SplitPolygon(Polygon3 *poly, const BspTree &tree) const;
+                                           const Plane3 &splitPlane,
+                       const AxisAlignedBox3 &box,
+                                           const float epsilon) const;
+
+        Polygon3 *SplitPolygon(Polygon3 *poly, const float epsilon) const;
 
         PolygonContainer mPolys;
@@ -270 +273 @@
         BspNode *GetFront();
 
-        Plane3 *GetPlane();
-
+        /** Returns split plane.
+        */
+        Plane3 GetPlane() const;
+
+        /** Replace front or back child with new child.
+        */
         void ReplaceChildLink(BspNode *oldChild, BspNode *newChild);
+        /** Replace front and back child.
+        */
         void SetupChildLinks(BspNode *b, BspNode *f);
-
-        /** Splits polygons with respect to the split plane.
-                @param polys the polygons to be split. the polygons are consumed and
-                           distributed to the containers frontPolys, backPolys, coincident.
-                @param frontPolys returns the polygons in the front of the split plane
-                @param backPolys returns the polygons in the back of the split plane
-                @param coincident returns the polygons coincident to the split plane
-
-                @returns the number of splits   
-        */
-        int SplitPolygons(PolygonContainer &polys, 
-                                          PolygonContainer &frontPolys, 
-                                          PolygonContainer &backPolys, 
-                                          PolygonContainer &coincident);
 
         friend ostream &operator<<(ostream &s, const BspInterior &A)
@@ -298 +293 @@
         /// Splitting plane corresponding to this node
         Plane3 mPlane;
+
         /// back node
         BspNode *mBack;
@@ -328 +324 @@
         void SetViewCell(BspViewCell *viewCell);
 
-        /** Adds ray sample contributions to the PVS.
-                @param sampleContributions the number contributions of the samples
-                @param contributingSampels the number of contributing rays
-                
-        */
-        void AddToPvs(const BoundedRayContainer &rays, 
-                                  int &sampleContributions,
-                                  int &contributingSamples);
+        VssRayContainer mVssRays;
 
 protected:
@@ -516 +505 @@
         BspViewCell *GetRootCell() const;
 
-        /** Parses the environment and stores the global BSP tree parameters
-        */
-        static void ParseEnvironment();
-
+        /** Returns epsilon of this tree.
+        */
+        float GetEpsilon() const;
 
 protected:
@@ -722 +710 @@
         bool SelectAxisAlignedPlane(Plane3 &plane, const PolygonContainer &polys) const;
 
-        /** Bounds ray and returns minT and maxT.
-                @returns true if ray hits BSP tree bounding box
-        */
-        bool BoundRay(const Ray &ray, float &minT, float &maxT) const;
-
         /** Subdivides the rays into front and back rays according to the split plane.
                 
@@ -768 +751 @@
         */
         float AccumulatedRayLength(BoundedRayContainer &rays) const;
+
+        /** Splits polygons with respect to the split plane.
+                @param polys the polygons to be split. the polygons are consumed and
+                           distributed to the containers frontPolys, backPolys, coincident.
+                @param frontPolys returns the polygons in the front of the split plane
+                @param backPolys returns the polygons in the back of the split plane
+                @param coincident returns the polygons coincident to the split plane
+
+                @returns the number of splits   
+        */
+        int SplitPolygons(const Plane3 &plane,
+                                          PolygonContainer &polys, 
+                                          PolygonContainer &frontPolys, 
+                                          PolygonContainer &backPolys, 
+                                          PolygonContainer &coincident) const;
+
+        /** Adds ray sample contributions to the PVS.
+                @param sampleContributions the number contributions of the samples
+                @param contributingSampels the number of contributing rays
+                
+        */
+        void AddToPvs(BspLeaf *leaf,
+                                  const BoundedRayContainer &rays, 
+                                  int &sampleContributions,
+                                  int &contributingSamples);
 
         /// Pointer to the root of the tree.
@@ -852 +860 @@
         bool mPvsUseArea;
 
+        /// epsilon where two points are still considered equal
+        float mEpsilon;
+
 private:
         

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.cpp

@@ -68 +68 @@
 
 void ViewCellsManager::ComputeSampleContributions(const RayContainer &rays, 
-                                                                                                  const bool castRays,
                                                                                                   int &sampleContributions,
-                                                                                                  int &contributingSamples)
+                                                                                                  int &contributingSamples,
+                                                                                                  const bool castRays)
 {
         // view cells not yet constructed
@@ -168 +168 @@
 SimulationStatistics ViewCellsManager::SimulateRendering() const
 {
-        if (!ViewCellsConstructed())
+        if (!ViewCellsConstructed() || !mRenderSimulator)
                 return SimulationStatistics();
 
@@ -240 +240 @@
                 return 0;
 
-        Debug << "Constructing bsp view cells" << endl;
-
         int sampleContributions = 0;
         
@@ -247 +245 @@
 
         int limit = min (mConstructionSamples, (int)rays.size());
+
+        // construct view cells using the collected samples
+        Debug << "construcing bsp view cells from " 
+                  << limit << " samples " << endl;
 
     for (int i = 0; i < limit; ++ i)
@@ -332 +334 @@
         int pvsSize = 0;
 
-        BspViewCellsStatistics stat;
-        mBspTree->EvaluateViewCellsStats(stat);
-        Debug << "original view cell partition:\n" << stat << endl;
-                
-        if (1) // export view cells
+        BspViewCellsStatistics vcStats;
+        mBspTree->EvaluateViewCellsStats(vcStats);
+        Debug << "original view cell partition:\n" << vcStats << endl;
+                
+        if (1) // export view cells before merge
         {
                 cout << "exporting initial view cells (=leaves) ... ";
@@ -344 +346 @@
                 {
                         exporter->SetWireframe();
-                        exporter->ExportBspLeaves(*mBspTree, stat.maxPvs);
+                        exporter->ExportBspLeaves(*mBspTree, vcStats.maxPvs);
                         //exporter->ExportBspViewCellPartition(*mBspTree, 0);
                                 
@@ -360 +362 @@
                 }
                 cout << "finished" << endl;
-        }
+
+                //-- render simulation
+                cout << "\nevaluating bsp view cells render time before merge ... ";
+                 
+                const SimulationStatistics ss = SimulateRendering();
+                         
+                cout << " finished" << endl;
+                cout << ss << endl;
+                Debug << ss << endl;
+        }
+
 
         cout << "starting post processing using " << mPostProcessSamples << " samples ... ";
-                
+
         long startTime = GetTime();
-
 
         //-- merge or subdivide view cells
@@ -408 +419 @@
         cout << "merged " << merged << " view cells in "
                  << TimeDiff(startTime, GetTime()) *1e-3 << " secs" << endl;
+
+        // statistics after post process
+        vcStats.Reset();
+        mBspTree->EvaluateViewCellsStats(vcStats);
+        Debug << "post processed view cell partition:\n" << vcStats << endl;
 
         return merged;
@@ -956 +972 @@
 mVspKdTree(vspKdTree)
 {
-   mRenderSimulator = NULL; // TODO: new VspKdRenderSimulator(vspKdTree);
+   mRenderSimulator = NULL;//new VspKdRenderSimulator(vspKdTree);
    InitRenderSimulator();
 }
@@ -983 +999 @@
                 return 0;
 
-        //if (castRay) // TODO
+        //if (castRay)
+                //mVspKdTree->CastRay(ray);
 
         int sampleContributions = 0;
 
-/*
-    for (int i = 0; i < (int)rays.size(); ++ i)
-        {
-                Ray *ray = rays[i];
-
-                mBspTree->CastRay(*ray);
-                        
-                Intersectable *term = 
-                        !ray->intersections.empty() ? ray->intersections[0].mObject : NULL
-                        ;
-                sampleContributions += 
-                        AddSampleContributions(*ray, ray->sourceObject.mObject, term);
-                        
-        }
-*/
         return sampleContributions;
 }
@@ -1121 +1123 @@
 mVspBspTree(vspBspTree)
 {
-   mRenderSimulator = NULL;// TODO new VspBspRenderSimulator(bspTree);
-   //InitRenderSimulator();
+        mRenderSimulator = new VspBspRenderSimulator(vspBspTree);
+        InitRenderSimulator();
 }
 
@@ -1182 +1184 @@
                 for (int j = 0; j < (int)ray.bspIntersections.size(); ++ j)
                 {       
-                        VspBspLeaf *leaf = NULL;  //TODO ray.bspIntersections[j].mLeaf;
+                        BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
                 
                         // if ray not in unbounded space
@@ -1206 +1208 @@
                 for (int j = 1; j < ((int)ray.bspIntersections.size() - 1); ++ j) 
                 {
-                        VspBspLeaf *leaf = NULL;// TODO = ray.bspIntersections[j].mLeaf;
+                        BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
 
                         if (leaf->GetViewCell() != mVspBspTree->GetRootCell())
@@ -1219 +1221 @@
                                                                                 const RayContainer &rays)
 {
-        return 0;// TODO
         if (!ViewCellsConstructed())
         {
@@ -1231 +1232 @@
         int pvsSize = 0;
 
-        VspBspViewCellsStatistics vcStats;
+        BspViewCellsStatistics vcStats;
         mVspBspTree->EvaluateViewCellsStats(vcStats);
         Debug << "original view cell partition:\n" << vcStats << endl;
@@ -1243 +1244 @@
                 {
                         exporter->SetWireframe();
-                        //exporter->ExportBspLeaves(*mVspBspTree, stat.maxPvs);
+                        exporter->ExportBspViewCellPartition(*mVspBspTree, vcStats.maxPvs);
                                 
                         if (0)
@@ -1284 +1285 @@
                 iit = ray->bspIntersections.begin();
 
-                VspBspLeaf *previousLeaf = NULL;//TODO(*iit).mLeaf;
+                BspLeaf *previousLeaf = (*iit).mLeaf;
                 ++ iit;
                 
                 for (; iit != ray->bspIntersections.end(); ++ iit)
                 {
-                        VspBspLeaf *leaf = NULL;//TODO(*iit).mLeaf;
+                        BspLeaf *leaf = (*iit).mLeaf;
 
                         if (ShouldMerge(leaf, previousLeaf))
@@ -1307 +1308 @@
                  << TimeDiff(startTime, GetTime()) *1e-3 << " secs" << endl;
 
+        // statistics after post process
+        vcStats.Reset();
+        mVspBspTree->EvaluateViewCellsStats(vcStats);
+        Debug << "post processed view cell partition:\n" << vcStats << endl;
+
         return merged;
 }
@@ -1319 +1325 @@
                                                                            const RayContainer &sampleRays)
 {
-        return; // TODO
         if (!ViewCellsConstructed())
                 return;
 
         //-- recount pvs
-        VspBspViewCellsStatistics vcStats;
+        BspViewCellsStatistics vcStats;
         mVspBspTree->EvaluateViewCellsStats(vcStats);
 
@@ -1330 +1335 @@
         {
                 cout << "exporting view cells after merge ... ";
-                /*Exporter *exporter = Exporter::GetExporter("merged_view_cells.x3d");
+                Exporter *exporter = Exporter::GetExporter("merged_view_cells.x3d");
 
                 if (exporter)
@@ -1337 +1342 @@
                         delete exporter;
                 }
-                */
+                
                 cout << "finished" << endl;
         }       
@@ -1343 +1348 @@
         //-- visualization of the BSP splits
         bool exportSplits = false;
-                environment->GetBoolValue("BspTree.Visualization.exportSplits", exportSplits);
+        environment->GetBoolValue("BspTree.Visualization.exportSplits", exportSplits);
         
         if (exportSplits)
@@ -1352 +1357 @@
         }
 
-//TODO  ExportVspBspPvs(objects, sampleRays);
+        ExportBspPvs(objects, sampleRays);
 }
 
@@ -1367 +1372 @@
                 exporter->SetForcedMaterial(m);
                 exporter->SetWireframe();
-                // TODO
-                //exporter->ExportBspSplits(*mVspBspTree, true);
+                
+                exporter->ExportBspSplits(*mVspBspTree, true);
 
                 // take forced material, else big scenes cannot be viewed
@@ -1412 +1417 @@
         const int raysOut = min((int)sampleRays.size(), mVisualizationSamples);
         cout << "visualization using " << mVisualizationSamples << " samples" << endl;
-        vector<Ray *> vcRays[leafOut];
-
-        if (0)
+        
+        if (1)
         {
                 //-- some random view cells and rays for output
-                vector<VspBspLeaf *> vspBspLeaves;
+                vector<BspLeaf *> vspBspLeaves;
 
                 for (int i = 0; i < leafOut; ++ i)
@@ -1424 +1428 @@
                 for (int i = 0; i < (int)vspBspLeaves.size(); ++ i)
                 {
+                        RayContainer vcRays;
                         cout << "creating output for view cell " << i << " ... ";
                         // check whether we can add the current ray to the output rays
-                        for (int k = 0; k < raysOut; ++ k) 
+                        /*for (int k = 0; k < raysOut; ++ k) 
                         {
                                 Ray *ray = sampleRays[k];
@@ -1436 +1441 @@
                                         if (vspBspLeaves[i]->GetViewCell() == leaf->GetViewCell()) 
                                         {
-                                                vcRays[i].push_back(ray);
+                                                vcRays.push_back(ray);
                                         }
                                 }
                         }
+                        */
 
                         Intersectable::NewMail();
@@ -1474 +1480 @@
 
                         Debug << i << ": pvs size=" << (int)vc->GetPvs().GetSize() 
-                                        << ", piercing rays=" << (int)vcRays[i].size() << endl;
+                                        << ", piercing rays=" << (int)vcRays.size() << endl;
 
                         // export rays piercing this view cell
-                        exporter->ExportRays(vcRays[i], 1000, RgbColor(0, 1, 0));
-
-                        m.mDiffuseColor = RgbColor(1, 0, 0);
+                        exporter->ExportRays(vcRays, 1000, RgbColor(1, 0, 0));
+                        exporter->ExportRays(vspBspLeaves[i]->mVssRays, RgbColor(1, 1, 1));
+                        
+                        m.mDiffuseColor = RgbColor(0, 1, 1);
                         exporter->SetForcedMaterial(m);
 
-                        // exporter->SetWireframe();
-                        exporter->SetFilled();
+                        exporter->SetWireframe();
+                        //exporter->SetFilled();
 
                         // output PVS of view cell
@@ -1527 +1534 @@
                 {
                         cout << "creating output for view cell " << i << " ... ";
-                        
+                        RayContainer vcRays;
             Intersectable::NewMail();
                         BspViewCell *vc = dynamic_cast<BspViewCell *>(viewCells[i]);
@@ -1539 +1546 @@
                                 for     (int j = 0; j < (int)ray->bspIntersections.size(); ++ j)
                                 {
-                                        VspBspLeaf *leaf = NULL; // TODO = ray->bspIntersections[j].mLeaf;
+                                        BspLeaf *leaf = ray->bspIntersections[j].mLeaf;
 
                                         if (vc == leaf->GetViewCell()) 
                                         {
-                                                vcRays[i].push_back(ray);
+                                                vcRays.push_back(ray);
                                         }
                                 }
@@ -1572 +1579 @@
                         
                         Debug << i << ": pvs size=" << (int)vc->GetPvs().GetSize() 
-                                        << ", piercing rays=" << (int)vcRays[i].size() << endl;
+                                        << ", piercing rays=" << (int)vcRays.size() << endl;
 
                         
                         // export rays piercing this view cell
-                        exporter->ExportRays(vcRays[i], 1000, RgbColor(0, 1, 0));
+                        exporter->ExportRays(vcRays, 1000, RgbColor(0, 1, 0));
         
                         m.mDiffuseColor = RgbColor(1, 0, 0);
@@ -1607 +1614 @@
 
 
-bool VspBspViewCellsManager::MergeVspBspLeafViewCells(VspBspLeaf *front, 
-                                                                                                          VspBspLeaf *back) const
+bool VspBspViewCellsManager::MergeVspBspLeafViewCells(BspLeaf *front, 
+                                                                                                          BspLeaf *back) const
 {
         BspViewCell *viewCell = 
@@ -1623 +1630 @@
         BspViewCell *bVc = back->GetViewCell();
 
-        // TODO
-        /*vector<VspBspLeaf *> fLeaves = fVc->mBspLeaves;
-        vector<VspBspLeaf *> bLeaves = bVc->mBspLeaves;
-
-        vector<VspBspLeaf *>::const_iterator it;
+        vector<BspLeaf *> fLeaves = fVc->mBspLeaves;
+        vector<BspLeaf *> bLeaves = bVc->mBspLeaves;
+
+        vector<BspLeaf *>::const_iterator it;
         
         for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
@@ -1642 +1648 @@
         DEL_PTR(fVc);
         DEL_PTR(bVc);
-*/
+
         return true;
 }
 
-bool VspBspViewCellsManager::ShouldMerge(VspBspLeaf *front, VspBspLeaf *back) const
+bool VspBspViewCellsManager::ShouldMerge(BspLeaf *front, BspLeaf *back) const
 {
         ViewCell *fvc = front->GetViewCell();

trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellsManager.h

@@ -20 +20 @@
 class VspKdTree;
 class AxisAlignedBox3;
-class VspBspLeaf;
+class BspLeaf;
 
 /**
@@ -57 +57 @@
         */
         void  ComputeSampleContributions(const RayContainer &rays, 
-                                                                         const bool castRays,
                                                                          int &sampleContributions,
-                                                                         int &contributingSamples);
+                                                                         int &contributingSamples,
+                                                                         const bool castRays = true);
 
 
@@ -350 +350 @@
         /** Merges view cells front and back leaf view cell.
         */
-        bool MergeVspBspLeafViewCells(VspBspLeaf *front, VspBspLeaf *back) const;
+        bool MergeVspBspLeafViewCells(BspLeaf *front, BspLeaf *back) const;
 
         /** Returns true if front and back leaf should be merged.
         */
-        bool ShouldMerge(VspBspLeaf *front, VspBspLeaf *back) const;
+        bool ShouldMerge(BspLeaf *front, BspLeaf *back) const;
 
         /// the view space partition BSP tree.

trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp

@@ -13 +13 @@
 #include "Exporter.h"
 #include "Plane3.h"
+#include "ViewCellBsp.h"
 
 //-- static members
@@ -24 +25 @@
 const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
 
-int VspBspNode::sMailId = 1;
 
 int VspBspTree::sFrontId = 0; 
@@ -30 +30 @@
 int VspBspTree::sFrontAndBackId = 0;
 
+
 /****************************************************************/
-/*                  class VspBspNode implementation                */
+/*                  class VspBspTree implementation             */
 /****************************************************************/
 
-VspBspNode::VspBspNode(): 
-mParent(NULL)
-{}
-
-VspBspNode::VspBspNode(VspBspInterior *parent): 
-mParent(parent)
-{}
-
-
-bool VspBspNode::IsRoot() const 
-{
-        return mParent == NULL;
-}
-
-VspBspInterior *VspBspNode::GetParent() 
-{ 
-        return mParent; 
-}
-
-void VspBspNode::SetParent(VspBspInterior *parent)
-{
-        mParent = parent;
-}
-
-/****************************************************************/
-/*              class VspBspInterior implementation                */
-/****************************************************************/
-
-
-VspBspInterior::VspBspInterior(const Plane3 &plane): 
-mPlane(plane), mFront(NULL), mBack(NULL)
-{}
-
-VspBspInterior::~VspBspInterior() 
-{ 
-        DEL_PTR(mFront); 
-        DEL_PTR(mBack);
-}
-
-bool VspBspInterior::IsLeaf() const
-{ 
-        return false; 
-}
-
-VspBspNode *VspBspInterior::GetBack() 
-{
-        return mBack;
-}
-
-VspBspNode *VspBspInterior::GetFront() 
-{
-        return mFront;
-}
-
-Plane3 *VspBspInterior::GetPlane()
-{
-        return &mPlane;
-}
-
-void VspBspInterior::ReplaceChildLink(VspBspNode *oldChild, VspBspNode *newChild) 
-{
-        if (mBack == oldChild)
-                mBack = newChild;
+VspBspTree::VspBspTree(): 
+mRoot(NULL),
+mPvsUseArea(true)
+{
+        mRootCell = new BspViewCell();
+
+        Randomize(); // initialise random generator for heuristics
+
+        //-- termination criteria for autopartition
+        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
+        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
+        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
+        environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);     
+        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
+        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
+
+        //-- factors for bsp tree split plane heuristics
+        environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
+        environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
+        environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
+
+        //-- termination criteria for axis aligned split
+        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAaCtDivCi);
+        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
+        
+        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays", 
+                                                         mTermMinRaysForAxisAligned);
+        
+        //-- partition criteria
+        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
+        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
+        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
+        environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mSplitBorder);
+
+        environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
+        environment->GetIntValue("VspBspTree.maxTests", mMaxTests);
+
+    Debug << "BSP max depth: " << mTermMaxDepth << endl;
+        Debug << "BSP min PVS: " << mTermMinPvs << endl;
+        Debug << "BSP min area: " << mTermMinArea << endl;
+        Debug << "BSP min rays: " << mTermMinRays << endl;
+        Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
+        Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;
+
+        Debug << "Split plane strategy: ";
+        if (mSplitPlaneStrategy & RANDOM_POLYGON)
+                Debug << "random polygon ";
+        if (mSplitPlaneStrategy & AXIS_ALIGNED)
+                Debug << "axis aligned ";
+        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+                Debug << "least ray splits ";
+        if (mSplitPlaneStrategy & BALANCED_RAYS)
+                Debug << "balanced rays ";
+        if (mSplitPlaneStrategy & PVS)
+                Debug << "pvs";
+
+        Debug << endl;
+}
+
+
+const BspTreeStatistics &VspBspTree::GetStatistics() const 
+{
+        return mStat;
+}
+
+
+VspBspTree::~VspBspTree()
+{
+        DEL_PTR(mRoot);
+        DEL_PTR(mRootCell);
+}
+
+int VspBspTree::AddMeshToPolygons(Mesh *mesh, 
+                                                                  PolygonContainer &polys, 
+                                                                  MeshInstance *parent)
+{
+        FaceContainer::const_iterator fi;
+        
+        // copy the face data to polygons
+        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
+        {
+                Polygon3 *poly = new Polygon3((*fi), mesh);
+                
+                if (poly->Valid(mEpsilon))
+                {
+                        poly->mParent = parent; // set parent intersectable
+                        polys.push_back(poly);
+                }
+                else
+                        DEL_PTR(poly);
+        }
+        return (int)mesh->mFaces.size();
+}
+
+int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells, 
+                                                          PolygonContainer &polys, 
+                                                          int maxObjects)
+{
+        int limit = (maxObjects > 0) ? 
+                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
+  
+        int polysSize = 0;
+
+        for (int i = 0; i < limit; ++ i)
+        {
+                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
+                {
+                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
+                        polysSize += 
+                                AddMeshToPolygons(viewCells[i]->GetMesh(), 
+                                                                  polys, 
+                                                                  viewCells[i]);
+                }
+        }
+
+        return polysSize;
+}
+
+int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects, 
+                                                                 PolygonContainer &polys, 
+                                                                 int maxObjects)
+{
+        int limit = (maxObjects > 0) ? 
+                Min((int)objects.size(), maxObjects) : (int)objects.size();
+  
+        for (int i = 0; i < limit; ++i)
+        {
+                Intersectable *object = objects[i];//*it;
+                Mesh *mesh = NULL;
+
+                switch (object->Type()) // extract the meshes
+                {
+                case Intersectable::MESH_INSTANCE:
+                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
+                        break;
+                case Intersectable::VIEW_CELL:
+                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
+                        break;
+                        // TODO: handle transformed mesh instances
+                default:
+                        Debug << "intersectable type not supported" << endl;
+                        break;
+                }
+                
+        if (mesh) // copy the mesh data to polygons
+                {
+                        mBox.Include(object->GetBox()); // add to BSP tree aabb
+                        AddMeshToPolygons(mesh, polys, mRootCell);
+                }
+        }
+
+        return (int)polys.size();
+}
+
+void VspBspTree::Construct(const VssRayContainer &sampleRays)
+{
+    mStat.nodes = 1;
+        mBox.Initialize();      // initialise BSP tree bounding box
+        
+        PolygonContainer polys;
+        RayInfoContainer *rays = new RayInfoContainer();
+
+        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
+
+        long startTime = GetTime();
+
+        Debug << "**** Extracting polygons from rays ****\n";
+
+        Intersectable::NewMail();
+
+        //-- extract polygons intersected by the rays
+        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
+        {
+                VssRay *ray = *rit;
+        
+                if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
+                {
+                        ray->mTerminationObject->Mail();
+                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
+                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
+                }
+
+                if (ray->mOriginObject && !ray->mOriginObject->Mailed())
+                {
+                        ray->mOriginObject->Mail();
+                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
+                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
+                }
+        }
+
+        // compute bounding box
+        Polygon3::IncludeInBox(polys, mBox);
+
+        //-- store rays
+        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
+        {
+                VssRay *ray = *rit;
+                
+                float minT, maxT;
+
+                // TODO: not very efficient to implictly cast between rays types ...
+                if (mBox.GetRaySegment(*ray, minT, maxT))
+                {
+                        float len = ray->Length();
+                        
+                        if (!len) 
+                                len = Limits::Small;
+                        
+                        rays->push_back(RayInfo(ray, minT / len, maxT / len));
+                }
+        }
+
+        mStat.polys = (int)polys.size();
+
+        Debug << "**** Finished polygon extraction ****" << endl;
+        Debug << (int)polys.size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
+        Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+
+        Construct(polys, rays);
+
+        // clean up polygons
+        CLEAR_CONTAINER(polys);
+}
+
+void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
+{
+        std::stack<VspBspTraversalData> tStack;
+
+        mRoot = new BspLeaf();
+
+        // constrruct root node geometry
+        BspNodeGeometry *geom = new BspNodeGeometry();
+        ConstructGeometry(mRoot, *geom);
+
+        VspBspTraversalData tData(mRoot, 
+                                                          new PolygonContainer(polys), 
+                                                          0,
+                                                          rays, 
+                              ComputePvsSize(*rays), 
+                                                          geom->GetArea(), 
+                                                          geom);
+
+        tStack.push(tData);
+
+        mStat.Start();
+        cout << "Contructing vsp bsp tree ... ";
+
+        long startTime = GetTime();
+        
+        while (!tStack.empty()) 
+        {
+                tData = tStack.top();
+
+            tStack.pop();
+
+                // subdivide leaf node
+                BspNode *r = Subdivide(tStack, tData);
+
+                if (r == mRoot)
+                        Debug << "BSP tree construction time spent at root: " 
+                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
+        }
+
+        cout << "finished\n";
+
+        mStat.Stop();
+}
+
+bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
+{
+        return 
+                (((int)data.mRays->size() <= mTermMinRays) ||
+                 (data.mPvs <= mTermMinPvs) ||
+                 (data.mArea <= mTermMinArea) ||
+                // (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
+                 (data.mDepth >= mTermMaxDepth));
+}
+
+BspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, 
+                                                           VspBspTraversalData &tData)
+{
+        //-- terminate traversal  
+        if (TerminationCriteriaMet(tData))              
+        {
+                BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
+        
+                BspViewCell *viewCell = new BspViewCell();
+                
+                leaf->SetViewCell(viewCell);
+                viewCell->mBspLeaves.push_back(leaf);
+
+                //-- add pvs
+                if (viewCell != mRootCell)
+                {
+                        int conSamp = 0, sampCon = 0;
+                        AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
+                        
+                        mStat.contributingSamples += conSamp;
+                        mStat.sampleContributions += sampCon;
+                }
+
+                EvaluateLeafStats(tData);
+        
+                //-- clean up
+                
+                DEL_PTR(tData.mPolygons);
+                DEL_PTR(tData.mRays);
+                DEL_PTR(tData.mGeometry);
+
+                return leaf;
+        }
+
+        //-- continue subdivision
+        PolygonContainer coincident;
+        
+        VspBspTraversalData tFrontData(new PolygonContainer(), 
+                                                                   tData.mDepth + 1,
+                                                                   new RayInfoContainer(), 
+                                                                   new BspNodeGeometry());
+
+        VspBspTraversalData tBackData(new PolygonContainer(), 
+                                                                  tData.mDepth + 1,
+                                                                  new RayInfoContainer(), 
+                                                                  new BspNodeGeometry());
+
+        // create new interior node and two leaf nodes
+        BspInterior *interior = 
+                SubdivideNode(tData, tFrontData, tBackData, coincident);
+
+#ifdef _DEBUG   
+        if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
+        {       for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
+                        Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
+                Debug << endl;}
+#endif
+
+        // push the children on the stack
+        tStack.push(tFrontData);
+        tStack.push(tBackData);
+
+        // cleanup
+        DEL_PTR(tData.mNode);   
+
+        DEL_PTR(tData.mPolygons);
+        DEL_PTR(tData.mRays);
+        DEL_PTR(tData.mGeometry);
+
+        return interior;
+}
+
+BspInterior *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
+                                                                           VspBspTraversalData &frontData,
+                                                                           VspBspTraversalData &backData,
+                                                                           PolygonContainer &coincident)
+{
+        mStat.nodes += 2;
+        
+        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
+        
+        // select subdivision plane
+        BspInterior *interior = new BspInterior(SelectPlane(leaf, tData)); 
+
+#ifdef _DEBUG
+        Debug << interior << endl;
+#endif
+        
+        // subdivide rays into front and back rays
+        SplitRays(interior->GetPlane(), 
+                          *tData.mRays, 
+                          *frontData.mRays, 
+                          *backData.mRays);
+        
+        // subdivide polygons with plane
+        mStat.splits += SplitPolygons(interior->GetPlane(),
+                                                                  *tData.mPolygons, 
+                                          *frontData.mPolygons, 
+                                                                  *backData.mPolygons,
+                                                                  coincident);
+
+    // compute pvs
+        frontData.mPvs = ComputePvsSize(*frontData.mRays);
+        backData.mPvs = ComputePvsSize(*backData.mRays);
+
+        // split geometry and compute area
+        if (1)
+        {
+                tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
+                                                                           *backData.mGeometry, 
+                                                                           interior->GetPlane(),
+                                                                           mBox,
+                                                                           mEpsilon);
+        
+                
+                frontData.mArea = frontData.mGeometry->GetArea();
+                backData.mArea = backData.mGeometry->GetArea();
+        }
+
+        // compute accumulated ray length
+        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
+        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
+
+        //-- create front and back leaf
+
+        BspInterior *parent = leaf->GetParent();
+
+        // replace a link from node's parent
+        if (!leaf->IsRoot())
+        {
+                parent->ReplaceChildLink(leaf, interior);
+                interior->SetParent(parent);
+        }
+        else // new root
+        {
+                mRoot = interior;
+        }
+
+        // and setup child links
+        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
+        
+        frontData.mNode = interior->GetFront();
+        backData.mNode = interior->GetBack();
+        
+        //DEL_PTR(leaf);
+        return interior;
+}
+
+void VspBspTree::AddToPvs(BspLeaf *leaf,
+                                                  const RayInfoContainer &rays, 
+                                                  int &sampleContributions,
+                                                  int &contributingSamples)
+{
+        sampleContributions = 0;
+        contributingSamples = 0;
+
+    RayInfoContainer::const_iterator it, it_end = rays.end();
+
+        ViewCell *vc = leaf->GetViewCell();
+
+        // add contributions from samples to the PVS
+        for (it = rays.begin(); it != it_end; ++ it)
+        {
+                int contribution = 0;
+                VssRay *ray = (*it).mRay;
+                        
+                if (ray->mTerminationObject)
+                        contribution += vc->GetPvs().AddSample(ray->mTerminationObject);
+                
+                if (ray->mOriginObject)
+                        contribution += vc->GetPvs().AddSample(ray->mOriginObject);
+
+                if (contribution)
+                {
+                        sampleContributions += contribution;
+                        ++ contributingSamples;
+                }
+                //leaf->mVssRays.push_back(ray);
+        }
+}
+
+void VspBspTree::SortSplitCandidates(const PolygonContainer &polys, 
+                                                                         const int axis, 
+                                                                         vector<SortableEntry> &splitCandidates) const
+{
+  splitCandidates.clear();
+  
+  int requestedSize = 2 * (int)polys.size();
+  // creates a sorted split candidates array  
+  splitCandidates.reserve(requestedSize);
+  
+  PolygonContainer::const_iterator it, it_end = polys.end();
+
+  AxisAlignedBox3 box;
+
+  // insert all queries 
+  for(it = polys.begin(); it != it_end; ++ it)
+  {
+          box.Initialize();
+          box.Include(*(*it));
+          
+          splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
+      splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
+  }
+  
+  stable_sort(splitCandidates.begin(), splitCandidates.end());
+}
+
+
+float VspBspTree::BestCostRatio(const PolygonContainer &polys,
+                                                                const AxisAlignedBox3 &box,
+                                                                const int axis,
+                                                                float &position,
+                                                                int &objectsBack,
+                                int &objectsFront) const
+{
+        vector<SortableEntry> splitCandidates;
+
+        SortSplitCandidates(polys, axis, splitCandidates);
+        
+        // go through the lists, count the number of objects left and right
+        // and evaluate the following cost funcion:
+        // C = ct_div_ci  + (ol + or)/queries
+        
+        int objectsLeft = 0, objectsRight = (int)polys.size();
+        
+        float minBox = box.Min(axis);
+        float maxBox = box.Max(axis);
+        float boxArea = box.SurfaceArea();
+  
+        float minBand = minBox + mSplitBorder * (maxBox - minBox);
+        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
+        
+        float minSum = 1e20f;
+        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
+
+        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci) 
+        {
+                switch ((*ci).type) 
+                {
+                        case SortableEntry::POLY_MIN:
+                                ++ objectsLeft;
+                                break;
+                        case SortableEntry::POLY_MAX:
+                            -- objectsRight;
+                                break;
+                        default:
+                                break;
+                }
+                
+                if ((*ci).value > minBand && (*ci).value < maxBand) 
+                {
+                        AxisAlignedBox3 lbox = box;
+                        AxisAlignedBox3 rbox = box;
+                        lbox.SetMax(axis, (*ci).value);
+                        rbox.SetMin(axis, (*ci).value);
+
+                        float sum = objectsLeft * lbox.SurfaceArea() + 
+                                                objectsRight * rbox.SurfaceArea();
+      
+                        if (sum < minSum) 
+                        {
+                                minSum = sum;
+                                position = (*ci).value;
+
+                                objectsBack = objectsLeft;
+                                objectsFront = objectsRight;
+                        }
+                }
+        }
+  
+        float oldCost = (float)polys.size();
+        float newCost = mAaCtDivCi + minSum / boxArea;
+        float ratio = newCost / oldCost;
+
+
+#if 0
+  Debug << "====================" << endl;
+  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
+      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
+#endif
+  return ratio;
+}
+
+bool VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
+                                        const PolygonContainer &polys) const
+{
+        AxisAlignedBox3 box;
+        box.Initialize();
+        
+        // create bounding box of region
+        Polygon3::IncludeInBox(polys, box);
+        
+        int objectsBack = 0, objectsFront = 0;
+        int axis = 0;
+        float costRatio = MAX_FLOAT;
+        Vector3 position;
+
+        //-- area subdivision
+        for (int i = 0; i < 3; ++ i) 
+        {
+                float p = 0;
+                float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
+                
+                if (r < costRatio)
+                {
+                        costRatio = r;
+                        axis = i;
+                        position = p;
+                }
+        }
+        
+        if (costRatio >= mMaxCostRatio)
+                return false;
+
+        Vector3 norm(0,0,0); norm[axis] = 1.0f;
+        plane = Plane3(norm, position);
+
+        return true;
+}
+
+Plane3 VspBspTree::SelectPlane(BspLeaf *leaf, VspBspTraversalData &data)
+{
+        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
+                ((int)data.mRays->size() > mTermMinRaysForAxisAligned))
+        {
+                Plane3 plane;
+                if (SelectAxisAlignedPlane(plane, *data.mPolygons)) // TODO: should be rays
+                        return plane;
+        }
+
+        // simplest strategy: just take next polygon
+        if (mSplitPlaneStrategy & RANDOM_POLYGON)
+        {
+        if (!data.mPolygons->empty())
+                {
+                        const int randIdx = Random((int)data.mPolygons->size());
+                        Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
+                        return nextPoly->GetSupportingPlane();
+                }
+                else
+                {
+                        //-- choose plane on midpoint of a ray
+                        const int candidateIdx = Random((int)data.mRays->size());
+                                                                        
+                        const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
+                        const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();
+
+                        const Vector3 pt = (maxPt + minPt) * 0.5;
+
+                        const Vector3 normal = (*data.mRays)[candidateIdx].mRay->GetDir();
+                        
+                        return Plane3(normal, pt);
+                }
+
+                return Plane3();
+        }
+
+        // use heuristics to find appropriate plane
+        return SelectPlaneHeuristics(leaf, data); 
+}
+
+Plane3 VspBspTree::SelectPlaneHeuristics(BspLeaf *leaf, VspBspTraversalData &data)
+{
+        float lowestCost = MAX_FLOAT;
+        Plane3 bestPlane;
+        Plane3 plane;
+
+        int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
+        
+        int candidateIdx = limit;
+        
+        for (int i = 0; i < limit; ++ i)
+        {
+                candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
+                
+                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
+
+                // evaluate current candidate
+                const float candidateCost = 
+                        SplitPlaneCost(poly->GetSupportingPlane(), data);
+
+                if (candidateCost < lowestCost)
+                {
+                        bestPlane = poly->GetSupportingPlane();
+                        lowestCost = candidateCost;
+                }
+        }
+        
+        //Debug << "lowest: " << lowestCost << endl;
+
+        //-- choose candidate planes extracted from rays
+        // we currently use different two methods chosen with
+        // equal probability
+        
+        // take 3 ray endpoints, where two are minimum and one a maximum
+        // point or the other way round
+        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
+        {
+                candidateIdx = Random((int)data.mRays->size());
+        
+                RayInfo rayInf = (*data.mRays)[candidateIdx];
+
+                const Vector3 minPt = rayInf.ExtrapOrigin();
+                const Vector3 maxPt = rayInf.ExtrapTermination();
+
+                const Vector3 pt = (maxPt + minPt) * 0.5;
+                const Vector3 normal = Normalize(rayInf.mRay->GetDir());
+
+                plane = Plane3(normal, pt);
+
+                const float candidateCost = SplitPlaneCost(plane, data);
+
+                if (candidateCost < lowestCost)
+                {
+                        bestPlane = plane;      
+                        lowestCost = candidateCost;
+                }
+        }
+
+        // take plane normal as plane normal and the midpoint of the ray.
+        // PROBLEM: does not resemble any point where visibility is likely to change
+        //Debug << "lowest2: " << lowestCost << endl;
+        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
+        {
+                Vector3 pt[3];
+                int idx[3];
+                int cmaxT = 0;
+                int cminT = 0;
+                bool chooseMin = false;
+
+                for (int j = 0; j < 3; j ++)
+                {
+                        idx[j] = Random((int)data.mRays->size() * 2);
+                                
+                        if (idx[j] >= (int)data.mRays->size())
+                        {
+                                idx[j] -= (int)data.mRays->size();
+                                
+                                chooseMin = (cminT < 2);
+                        }
+                        else
+                                chooseMin = (cmaxT < 2);
+
+                        RayInfo rayInf = (*data.mRays)[idx[j]];
+                        pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
+                }       
+                        
+                plane = Plane3(pt[0], pt[1], pt[2]);
+
+                const float candidateCost = SplitPlaneCost(plane, data);
+
+                if (candidateCost < lowestCost)
+                {
+                        //Debug << "choose ray plane 2: " << candidateCost << endl;
+                        bestPlane = plane;
+                        
+                        lowestCost = candidateCost;
+                }
+        }       
+
+#ifdef _DEBUG
+        Debug << "plane lowest cost: " << lowestCost << endl;
+#endif
+        return bestPlane;
+}
+
+int VspBspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
+{
+        const int candidateIdx = Random(currentIdx --);
+
+        // swap candidates to avoid testing same plane 2 times
+        std::swap(polys[currentIdx], polys[candidateIdx]);
+        
+        return currentIdx;
+        //return Random((int)polys.size());
+}
+
+
+inline void VspBspTree::GenerateUniqueIdsForPvs()
+{
+        Intersectable::NewMail(); sBackId = ViewCell::sMailId;
+        Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
+        Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
+}
+
+float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
+                                                                 const VspBspTraversalData &data)
+{
+        float val = 0;
+
+        float sumBalancedRays = 0;
+        float sumRaySplits = 0;
+        
+        int frontPvs = 0;
+        int backPvs = 0;
+
+        // probability that view point lies in child
+        float pOverall = 0;
+        float pFront = 0;
+        float pBack = 0;
+
+        const bool pvsUseLen = false;
+
+        if (mSplitPlaneStrategy & PVS)
+        {
+                // create unique ids for pvs heuristics
+                GenerateUniqueIdsForPvs();
+        
+                if (mPvsUseArea) // use front and back cell areas to approximate volume
+                {       
+                        // construct child geometry with regard to the candidate split plane
+                        BspNodeGeometry frontCell;
+                        BspNodeGeometry backCell;
+                
+                        data.mGeometry->SplitGeometry(frontCell, 
+                                                                                  backCell, 
+                                                                                  candidatePlane,
+                                                                                  mBox,
+                                                                                  mEpsilon);
+                
+                        pFront = frontCell.GetArea();
+                        pBack = backCell.GetArea();
+
+                        pOverall = data.mArea;
+                }
+        }
+                
+        int limit;
+        bool useRand;
+
+        // choose test polyongs randomly if over threshold
+        if ((int)data.mRays->size() > mMaxTests)
+        {
+                useRand = true;
+                limit = mMaxTests;
+        }
         else
-                mFront = newChild;
-}
-
-void VspBspInterior::SetupChildLinks(VspBspNode *b, VspBspNode *f) 
-{
-    mBack = b;
-    mFront = f;
-}
-
-int VspBspInterior::SplitPolygons(PolygonContainer &polys, 
-                                                                  PolygonContainer &frontPolys, 
-                                                                  PolygonContainer &backPolys, 
-                                                                  PolygonContainer &coincident)
+        {
+                useRand = false;
+                limit = (int)data.mRays->size();
+        }
+
+        for (int i = 0; i < limit; ++ i)
+        {
+                const int testIdx = useRand ? Random(limit) : i;
+                RayInfo rayInf = (*data.mRays)[testIdx];
+
+                VssRay *ray = rayInf.mRay;
+                const int cf = rayInf.ComputeRayIntersection(candidatePlane, ray->mT);
+
+                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+                {
+                        sumBalancedRays += cf;
+                }
+                
+                if (mSplitPlaneStrategy & BALANCED_RAYS)
+                {
+                        if (cf == 0)
+                                ++ sumRaySplits;
+                }
+
+                if (mSplitPlaneStrategy & PVS)
+                {
+                        // in case the ray intersects an object
+                        // assure that we only count the object 
+                        // once for the front and once for the back side of the plane
+                        
+                        // add the termination object
+                        AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
+                        
+                        // add the source object
+                        AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
+                        
+                        // use number or length of rays to approximate volume
+                        if (!mPvsUseArea)
+                        {
+                                float len = 1;
+
+                                if (pvsUseLen) // use length of rays
+                                        len = rayInf.SqrSegmentLength();
+                        
+                                pOverall += len;
+
+                                if (cf == 1)
+                                        pFront += len;
+                                if (cf == -1)
+                                        pBack += len;
+                                if (cf == 0)
+                                {
+                                        // use length of rays to approximate volume
+                                        if (pvsUseLen) 
+                                        {
+                                                float newLen = len * 
+                                                        (rayInf.GetMaxT() - rayInf.mRay->mT) / 
+                                                        (rayInf.GetMaxT() - rayInf.GetMinT());
+                
+                                                if (candidatePlane.Side(rayInf.ExtrapOrigin()) <= 0)
+                                                {
+                                                        pBack += newLen;
+                                                        pFront += len - newLen;
+                                                }
+                                                else
+                                                {
+                                                        pFront += newLen;
+                                                        pBack += len - newLen;
+                                                }
+                                        }
+                                        else
+                                        {
+                                                ++ pFront;
+                                                ++ pBack;
+                                        }
+                                }
+                        }
+                }
+        }
+
+        const float raysSize = (float)data.mRays->size() + Limits::Small;
+
+        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
+                val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
+
+        if (mSplitPlaneStrategy & BALANCED_RAYS)
+                val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
+
+        const float denom = pOverall * (float)data.mPvs * 2.0f + Limits::Small;
+
+        if (mSplitPlaneStrategy & PVS)
+        {
+                val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
+
+                // give penalty to unbalanced split
+                if (1)
+                if (((pFront * 0.2 + Limits::Small) > pBack) || 
+                        (pFront < (pBack * 0.2 + Limits::Small)))
+                        val += 0.5;
+        }
+
+#ifdef _DEBUG
+        Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
+                  << " frontpvs: " << frontPvs << " pFront: " << pFront 
+                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
+#endif
+        return val;
+}
+
+void VspBspTree::AddObjToPvs(Intersectable *obj,
+                                                         const int cf,
+                                                         int &frontPvs,
+                                                         int &backPvs) const
+{
+        if (!obj)
+                return;
+        // TODO: does this really belong to no pvs?
+        //if (cf == Ray::COINCIDENT) return;
+
+        // object belongs to both PVS
+        if (cf >= 0)
+        {
+                if ((obj->mMailbox != sFrontId) && 
+                        (obj->mMailbox != sFrontAndBackId))
+                {
+                        ++ frontPvs;
+
+                        if (obj->mMailbox == sBackId)
+                                obj->mMailbox = sFrontAndBackId;        
+                        else
+                                obj->mMailbox = sFrontId;                                                               
+                }
+        }
+        
+        if (cf <= 0)
+        {
+                if ((obj->mMailbox != sBackId) &&
+                        (obj->mMailbox != sFrontAndBackId))
+                {
+                        ++ backPvs;
+
+                        if (obj->mMailbox == sFrontId)
+                                obj->mMailbox = sFrontAndBackId; 
+                        else
+                                obj->mMailbox = sBackId;                                
+                }
+        }
+}
+
+void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
+{
+        stack<BspNode *> nodeStack;
+        nodeStack.push(mRoot);
+ 
+        while (!nodeStack.empty()) 
+        {
+                BspNode *node = nodeStack.top();
+    
+                nodeStack.pop();
+    
+                if (node->IsLeaf()) 
+                {
+                        BspLeaf *leaf = (BspLeaf *)node;                
+                        leaves.push_back(leaf);
+                } 
+                else 
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+                        nodeStack.push(interior->GetBack());
+                        nodeStack.push(interior->GetFront());
+                }
+        }
+}
+
+AxisAlignedBox3 VspBspTree::GetBoundingBox() const
+{
+        return mBox;
+}
+
+BspNode *VspBspTree::GetRoot() const
+{
+        return mRoot;
+}
+
+void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
+{
+        // the node became a leaf -> evaluate stats for leafs
+        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
+
+        // store maximal and minimal depth
+        if (data.mDepth > mStat.maxDepth)
+                mStat.maxDepth = data.mDepth;
+
+        if (data.mDepth < mStat.minDepth)
+                mStat.minDepth = data.mDepth;
+
+        // accumulate depth to compute average depth
+        mStat.accumDepth += data.mDepth;
+        
+
+        if (data.mDepth >= mTermMaxDepth)
+                ++ mStat.maxDepthNodes;
+
+        if (data.mPvs < mTermMinPvs)
+                ++ mStat.minPvsNodes;
+
+        if ((int)data.mRays->size() < mTermMinRays)
+                ++ mStat.minRaysNodes;
+
+        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
+                ++ mStat.maxRayContribNodes;
+        
+        if (data.mGeometry->GetArea() <= mTermMinArea) 
+                ++ mStat.minAreaNodes;
+
+#ifdef _DEBUG
+        Debug << "BSP stats: "
+                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
+                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
+                  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
+                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), " 
+                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
+                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
+#endif
+}
+
+int VspBspTree::CastRay(Ray &ray)
+{
+        int hits = 0;
+  
+        stack<BspRayTraversalData> tStack;
+  
+        float maxt, mint;
+
+        if (!mBox.GetRaySegment(ray, mint, maxt))
+                return 0;
+
+        Intersectable::NewMail();
+
+        Vector3 entp = ray.Extrap(mint);
+        Vector3 extp = ray.Extrap(maxt);
+  
+        BspNode *node = mRoot;
+        BspNode *farChild = NULL;
+        
+        while (1)
+        {
+                if (!node->IsLeaf()) 
+                {
+                        BspInterior *in = dynamic_cast<BspInterior *>(node);
+
+                        Plane3 splitPlane = in->GetPlane();
+                        const int entSide = splitPlane.Side(entp);
+                        const int extSide = splitPlane.Side(extp);
+
+                        if (entSide < 0)
+                        {
+                                node = in->GetBack();
+
+                                if(extSide <= 0) // plane does not split ray => no far child
+                                        continue;
+                                        
+                                farChild = in->GetFront(); // plane splits ray
+
+                        } else if (entSide > 0)
+                        {
+                                node = in->GetFront();
+
+                                if (extSide >= 0) // plane does not split ray => no far child
+                                        continue;
+
+                                farChild = in->GetBack(); // plane splits ray                   
+                        }
+                        else // ray and plane are coincident
+                        {
+                                // WHAT TO DO IN THIS CASE ?
+                                //break;
+                                node = in->GetFront();
+                                continue;
+                        }
+
+                        // push data for far child
+                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
+
+                        // find intersection of ray segment with plane
+                        float t;
+                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
+                        maxt *= t;
+                        
+                } else // reached leaf => intersection with view cell
+                {
+                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
+      
+                        if (!leaf->GetViewCell()->Mailed())
+                        {
+                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
+                                leaf->GetViewCell()->Mail();
+                                ++ hits;
+                        }
+                        
+                        //-- fetch the next far child from the stack
+                        if (tStack.empty())
+                                break;
+      
+                        entp = extp;
+                        mint = maxt; // NOTE: need this?
+
+                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
+                                break;
+
+                        BspRayTraversalData &s = tStack.top();
+
+                        node = s.mNode;
+                        extp = s.mExitPoint;
+                        maxt = s.mMaxT;
+
+                        tStack.pop();
+                }
+        }
+
+        return hits;
+}
+
+bool VspBspTree::Export(const string filename)
+{
+        Exporter *exporter = Exporter::GetExporter(filename);
+
+        if (exporter) 
+        {
+                //exporter->ExportVspBspTree(*this);
+                return true;
+        }       
+
+        return false;
+}
+
+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 (!viewCell->Mailed()) 
+                        {
+                                viewCell->Mail();
+                                viewCells.push_back(viewCell);
+                        }
+                }
+                else 
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+                        nodeStack.push(interior->GetFront());
+                        nodeStack.push(interior->GetBack());
+                }
+        }
+}
+
+void VspBspTree::EvaluateViewCellsStats(BspViewCellsStatistics &stat) const
+{
+        stat.Reset();
+
+        stack<BspNode *> nodeStack;
+        nodeStack.push(mRoot);
+
+        ViewCell::NewMail();
+
+        // exclude root cell
+        mRootCell->Mail();
+
+        while (!nodeStack.empty()) 
+        {
+                BspNode *node = nodeStack.top();
+                nodeStack.pop();
+
+                if (node->IsLeaf()) 
+                {
+                        ++ stat.bspLeaves;
+
+                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
+
+                        if (!viewCell->Mailed()) 
+                        {
+                                viewCell->Mail();
+                                
+                                ++ stat.viewCells;
+                                const int pvsSize = viewCell->GetPvs().GetSize();
+
+                stat.pvs += pvsSize;
+
+                                if (pvsSize < 1)
+                                        ++ stat.emptyPvs;
+
+                                if (pvsSize > stat.maxPvs)
+                                        stat.maxPvs = pvsSize;
+
+                                if (pvsSize < stat.minPvs)
+                                        stat.minPvs = pvsSize;
+
+                                if ((int)viewCell->mBspLeaves.size() > stat.maxBspLeaves)
+                                        stat.maxBspLeaves = (int)viewCell->mBspLeaves.size();           
+                        }
+                }
+                else 
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+                        nodeStack.push(interior->GetFront());
+                        nodeStack.push(interior->GetBack());
+                }
+        }
+}
+
+BspTreeStatistics &VspBspTree::GetStat()
+{
+        return mStat;
+}
+
+float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
+{
+        float len = 0;
+
+        RayInfoContainer::const_iterator it, it_end = rays.end();
+
+        for (it = rays.begin(); it != it_end; ++ it)
+                len += (*it).SegmentLength();
+
+        return len;
+}
+
+int VspBspTree::SplitRays(const Plane3 &plane,
+                                                  RayInfoContainer &rays, 
+                                                  RayInfoContainer &frontRays, 
+                                                  RayInfoContainer &backRays)
+{
+        int splits = 0;
+
+        while (!rays.empty())
+        {
+                RayInfo bRay = rays.back();
+                VssRay *ray = bRay.mRay;
+
+                rays.pop_back();
+        
+                // get classification and receive new t
+                const int cf = bRay.ComputeRayIntersection(plane, ray->mT);
+
+                switch (cf)
+                {
+                case -1:
+                        backRays.push_back(bRay);
+                        break;
+                case 1:
+                        frontRays.push_back(bRay);
+                        break;
+                case 0: 
+                        //-- split ray
+
+                        // if start point behind plane
+                        if (plane.Side(bRay.ExtrapOrigin()) <= 0)
+                        {       
+                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
+                                frontRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
+                        }
+                        else
+                        {
+                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
+                                backRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
+                        }
+                        break;
+                default:
+                        Debug << "Should not come here 4" << endl;
+                        break;
+                }
+        }
+
+        return splits;
+}
+
+void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
+{
+        BspNode *lastNode;
+
+        do
+        {
+                lastNode = n;
+
+                // want to get planes defining geometry of this node => don't take
+                // split plane of node itself
+                n = n->GetParent();
+                
+                if (n)
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
+                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
+
+            if (interior->GetFront() != lastNode)
+                                halfSpace.ReverseOrientation();
+
+                        halfSpaces.push_back(halfSpace);
+                }
+        }
+        while (n);
+}
+
+void VspBspTree::ConstructGeometry(BspNode *n, 
+                                                                   BspNodeGeometry &cell) const
+{
+        PolygonContainer polys;
+        ConstructGeometry(n, polys);
+        cell.mPolys = polys;
+}
+
+void VspBspTree::ConstructGeometry(BspNode *n, 
+                                                                   PolygonContainer &cell) const
+{
+        vector<Plane3> halfSpaces;
+        ExtractHalfSpaces(n, halfSpaces);
+
+        PolygonContainer candidatePolys;
+
+        // bounded planes are added to the polygons (reverse polygons
+        // as they have to be outfacing
+        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
+        {
+                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
+                
+                if (p->Valid(mEpsilon))
+                {
+                        candidatePolys.push_back(p->CreateReversePolygon());
+                        DEL_PTR(p);
+                }
+        }
+
+        // add faces of bounding box (also could be faces of the cell)
+        for (int i = 0; i < 6; ++ i)
+        {
+                VertexContainer vertices;
+        
+                for (int j = 0; j < 4; ++ j)
+                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
+
+                candidatePolys.push_back(new Polygon3(vertices));
+        }
+
+        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
+        {
+                // polygon is split by all other planes
+                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
+                {
+                        if (i == j) // polygon and plane are coincident
+                                continue;
+
+                        VertexContainer splitPts;
+                        Polygon3 *frontPoly, *backPoly;
+
+                        const int cf = 
+                                candidatePolys[i]->ClassifyPlane(halfSpaces[j],
+                                                                                                 mEpsilon);
+                        
+                        switch (cf)
+                        {
+                                case Polygon3::SPLIT:
+                                        frontPoly = new Polygon3();
+                                        backPoly = new Polygon3();
+
+                                        candidatePolys[i]->Split(halfSpaces[j], 
+                                                                                         *frontPoly, 
+                                                                                         *backPoly,
+                                                                                         mEpsilon);
+
+                                        DEL_PTR(candidatePolys[i]);
+
+                                        if (frontPoly->Valid(mEpsilon))
+                                                candidatePolys[i] = frontPoly;
+                                        else
+                                                DEL_PTR(frontPoly);
+
+                                        DEL_PTR(backPoly);
+                                        break;
+                                case Polygon3::BACK_SIDE:
+                                        DEL_PTR(candidatePolys[i]);
+                                        break;
+                                // just take polygon as it is
+                                case Polygon3::FRONT_SIDE:
+                                case Polygon3::COINCIDENT:
+                                default:
+                                        break;
+                        }
+                }
+                
+                if (candidatePolys[i])
+                        cell.push_back(candidatePolys[i]);
+        }
+}
+
+void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &vcGeom) const
+{
+        vector<BspLeaf *> leaves = vc->mBspLeaves;
+
+        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
+
+        for (it = leaves.begin(); it != it_end; ++ it)
+                ConstructGeometry(*it, vcGeom);
+}
+
+int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors, 
+                                                   const bool onlyUnmailed) const
+{
+        PolygonContainer cell;
+
+        ConstructGeometry(n, cell);
+
+        stack<BspNode *> nodeStack;
+        nodeStack.push(mRoot);
+                
+        // planes needed to verify that we found neighbor leaf.
+        vector<Plane3> halfSpaces;
+        ExtractHalfSpaces(n, halfSpaces);
+
+        while (!nodeStack.empty()) 
+        {
+                BspNode *node = nodeStack.top();
+                nodeStack.pop();
+
+                if (node->IsLeaf())
+                {
+            if (node != n && (!onlyUnmailed || !node->Mailed()))
+                        {
+                                // test all planes of current node if candidate really
+                                // is neighbour
+                                PolygonContainer neighborCandidate;
+                                ConstructGeometry(node, neighborCandidate);
+                                
+                                bool isAdjacent = true;
+                                for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
+                                {
+                                        const int cf = 
+                                                Polygon3::ClassifyPlane(neighborCandidate, 
+                                                                                                halfSpaces[i],
+                                                                                                mEpsilon);
+
+                                        if (cf == Polygon3::BACK_SIDE)
+                                                isAdjacent = false;
+                                }
+
+                                if (isAdjacent)
+                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
+
+                                CLEAR_CONTAINER(neighborCandidate);
+                        }
+                }
+                else 
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
+        
+                        const int cf = Polygon3::ClassifyPlane(cell, 
+                                                                                                   interior->GetPlane(), 
+                                                                                                   mEpsilon);
+
+                        if (cf == Polygon3::FRONT_SIDE)
+                                nodeStack.push(interior->GetFront());
+                        else
+                                if (cf == Polygon3::BACK_SIDE)
+                                        nodeStack.push(interior->GetBack());
+                                else 
+                                {
+                                        // random decision
+                                        nodeStack.push(interior->GetBack());
+                                        nodeStack.push(interior->GetFront());
+                                }
+                }
+        }
+        
+        CLEAR_CONTAINER(cell);
+        return (int)neighbors.size();
+}
+
+BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
+{
+    stack<BspNode *> nodeStack;
+        nodeStack.push(mRoot);
+        
+        int mask = rand();
+  
+        while (!nodeStack.empty()) 
+        {
+                BspNode *node = nodeStack.top();
+                nodeStack.pop();
+          
+                if (node->IsLeaf()) 
+                {
+                        return dynamic_cast<BspLeaf *>(node);
+                } 
+                else 
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
+                        
+                        BspNode *next;
+        
+                        PolygonContainer cell;
+
+                        // todo: not very efficient: constructs full cell everytime
+                        ConstructGeometry(interior, cell);
+
+                        const int cf = Polygon3::ClassifyPlane(cell, halfspace, mEpsilon);
+
+                        if (cf == Polygon3::BACK_SIDE)
+                                next = interior->GetFront();
+                        else
+                                if (cf == Polygon3::FRONT_SIDE)
+                                        next = interior->GetFront();
+                        else 
+                        {
+                                // random decision
+                                if (mask & 1)
+                                        next = interior->GetBack();
+                                else
+                                        next = interior->GetFront();
+                                mask = mask >> 1;
+                        }
+
+                        nodeStack.push(next);
+                }
+        }
+        
+        return NULL;
+}
+
+BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
+{
+        stack<BspNode *> nodeStack;
+        
+        nodeStack.push(mRoot);
+
+        int mask = rand();
+        
+        while (!nodeStack.empty()) 
+        {
+                BspNode *node = nodeStack.top();
+                nodeStack.pop();
+                
+                if (node->IsLeaf()) 
+                {
+                        if ( (!onlyUnmailed || !node->Mailed()) )
+                                return dynamic_cast<BspLeaf *>(node);
+                }
+                else 
+                {
+                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
+
+                        // random decision
+                        if (mask & 1)
+                                nodeStack.push(interior->GetBack());
+                        else
+                                nodeStack.push(interior->GetFront());
+
+                        mask = mask >> 1;
+                }
+        }
+        
+        return NULL;
+}
+
+int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
+{
+        int pvsSize = 0;
+
+        RayInfoContainer::const_iterator rit, rit_end = rays.end();
+
+        Intersectable::NewMail();
+
+        for (rit = rays.begin(); rit != rays.end(); ++ rit)
+        {
+                VssRay *ray = (*rit).mRay;
+                
+                if (ray->mOriginObject)
+                {
+                        if (!ray->mOriginObject->Mailed())
+                        {
+                                ray->mOriginObject->Mail();
+                                ++ pvsSize;
+                        }
+                }
+                if (ray->mTerminationObject)
+                {
+                        if (!ray->mTerminationObject->Mailed())
+                        {
+                                ray->mTerminationObject->Mail();
+                                ++ pvsSize;
+                        }
+                }
+        }
+
+        return pvsSize;
+}
+
+float VspBspTree::GetEpsilon() const
+{
+        return mEpsilon;
+}
+
+BspViewCell *VspBspTree::GetRootCell() const
+{
+        return mRootCell;
+}
+
+int VspBspTree::SplitPolygons(const Plane3 &plane,
+                                                          PolygonContainer &polys, 
+                                                          PolygonContainer &frontPolys, 
+                                                          PolygonContainer &backPolys, 
+                                                          PolygonContainer &coincident) const
 {
         int splits = 0;
@@ -120 +1683 @@
 
                 // classify polygon
-                const int cf = poly->ClassifyPlane(mPlane);
+                const int cf = poly->ClassifyPlane(plane, mEpsilon);
 
                 switch (cf)
@@ -146 +1709 @@
         return splits;
 }
-
-/****************************************************************/
-/*                  class VspBspLeaf implementation                */
-/****************************************************************/
-VspBspLeaf::VspBspLeaf(): mViewCell(NULL)
-{
-}
-
-VspBspLeaf::VspBspLeaf(BspViewCell *viewCell): 
-mViewCell(viewCell)
-{
-}
-
-VspBspLeaf::VspBspLeaf(VspBspInterior *parent): 
-VspBspNode(parent), mViewCell(NULL)
-{}
-
-VspBspLeaf::VspBspLeaf(VspBspInterior *parent, BspViewCell *viewCell): 
-VspBspNode(parent), mViewCell(viewCell)
-{
-}
-
-BspViewCell *VspBspLeaf::GetViewCell() const
-{
-        return mViewCell;
-}
-
-void VspBspLeaf::SetViewCell(BspViewCell *viewCell)
-{
-        mViewCell = viewCell;
-}
-
-bool VspBspLeaf::IsLeaf() const 
-{ 
-        return true; 
-}
-
-void VspBspLeaf::AddToPvs(const RayInfoContainer &rays, 
-                                                  int &sampleContributions,
-                                                  int &contributingSamples,
-                                                  bool storeLeavesWithRays)
-{
-        sampleContributions = 0;
-        contributingSamples = 0;
-
-    RayInfoContainer::const_iterator it, it_end = rays.end();
-
-        // add contributions from samples to the PVS
-        for (it = rays.begin(); it != it_end; ++ it)
-        {
-                int contribution = 0;
-                VssRay *ray = (*it).mRay;
-                        
-                if (ray->mTerminationObject)
-                        contribution += mViewCell->GetPvs().AddSample(ray->mTerminationObject);
-                
-                if (ray->mOriginObject)
-                        contribution += mViewCell->GetPvs().AddSample(ray->mOriginObject);
-
-                if (contribution > 0)
-                {
-                        sampleContributions += contribution;
-                        ++ contributingSamples;
-                }
-        }
-}
-
-
-/****************************************************************/
-/*                  class VspBspTree implementation             */
-/****************************************************************/
-
-VspBspTree::VspBspTree(): 
-mRoot(NULL),
-mPvsUseArea(true)
-{
-        mRootCell = new BspViewCell();
-
-        Randomize(); // initialise random generator for heuristics
-
-        //-- termination criteria for autopartition
-        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
-        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
-        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
-        environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);     
-        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
-        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
-
-        //-- factors for bsp tree split plane heuristics
-        environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
-        environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
-        environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
-
-        //-- termination criteria for axis aligned split
-        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAaCtDivCi);
-        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
-        
-        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays", 
-                                                         mTermMinRaysForAxisAligned);
-        
-        //-- partition criteria
-        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
-        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
-        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
-        environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mSplitBorder);
-
-        environment->GetFloatValue("VspBspTree.Construction.sideTolerance", Vector3::sDistTolerance);
-        Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
-
-    Debug << "BSP max depth: " << mTermMaxDepth << endl;
-        Debug << "BSP min PVS: " << mTermMinPvs << endl;
-        Debug << "BSP min area: " << mTermMinArea << endl;
-        Debug << "BSP min rays: " << mTermMinRays << endl;
-        Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
-        Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;
-
-        Debug << "Split plane strategy: ";
-        if (mSplitPlaneStrategy & RANDOM_POLYGON)
-                Debug << "random polygon ";
-        if (mSplitPlaneStrategy & AXIS_ALIGNED)
-                Debug << "axis aligned ";
-        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-                Debug << "least ray splits ";
-        if (mSplitPlaneStrategy & BALANCED_RAYS)
-                Debug << "balanced rays ";
-        if (mSplitPlaneStrategy & PVS)
-                Debug << "pvs";
-
-        Debug << endl;
-}
-
-
-const VspBspTreeStatistics &VspBspTree::GetStatistics() const 
-{
-        return mStat;
-}
-
-void VspBspTreeStatistics::Print(ostream &app) const
-{
-        app << "===== VspBspTree statistics ===============\n";
-
-        app << setprecision(4);
-
-        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
-
-        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
-
-        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
-
-        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
-
-        app << "#N_SPLITS ( Number of splits )\n" << splits << "\n";
-
-        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
-                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
-
-        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
-                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
-
-        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n" 
-                << minPvsNodes * 100 / (double)Leaves() << endl;
-
-        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n" 
-                <<      minRaysNodes * 100 / (double)Leaves() << endl;
-
-        app << "#N_PMINAREALEAVES  ( Percentage of leaves with mininum area )\n"
-                << minAreaNodes * 100 / (double)Leaves() << endl;
-
-        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n" 
-                <<      maxRayContribNodes * 100 / (double)Leaves() << endl;
-
-        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
-
-        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
-
-        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
-
-        app << "#N_INPUT_POLYGONS (number of input polygons )\n" <<     polys << endl;
-
-        //app << "#N_PVS: " << pvs << endl;
-
-        app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
-                 (polys + splits) / (double)polys << endl;
-        
-        app << "===== END OF VspBspTree statistics ==========\n";
-}
-
-
-VspBspTree::~VspBspTree()
-{
-        DEL_PTR(mRoot);
-        DEL_PTR(mRootCell);
-}
-
-int VspBspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
-{
-        FaceContainer::const_iterator fi;
-        
-        // copy the face data to polygons
-        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
-        {
-                Polygon3 *poly = new Polygon3((*fi), mesh);
-                
-                if (poly->Valid())
-                {
-                        poly->mParent = parent; // set parent intersectable
-                        polys.push_back(poly);
-                }
-                else
-                        DEL_PTR(poly);
-        }
-        return (int)mesh->mFaces.size();
-}
-
-int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells, 
-                                                          PolygonContainer &polys, 
-                                                          int maxObjects)
-{
-        int limit = (maxObjects > 0) ? 
-                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
-  
-        int polysSize = 0;
-
-        for (int i = 0; i < limit; ++ i)
-        {
-                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
-                {
-                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
-                        polysSize += AddMeshToPolygons(viewCells[i]->GetMesh(), polys, viewCells[i]);
-                }
-        }
-
-        return polysSize;
-}
-
-int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects, PolygonContainer &polys, int maxObjects)
-{
-        int limit = (maxObjects > 0) ? Min((int)objects.size(), maxObjects) : (int)objects.size();
-  
-        for (int i = 0; i < limit; ++i)
-        {
-                Intersectable *object = objects[i];//*it;
-                Mesh *mesh = NULL;
-
-                switch (object->Type()) // extract the meshes
-                {
-                case Intersectable::MESH_INSTANCE:
-                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
-                        break;
-                case Intersectable::VIEW_CELL:
-                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
-                        break;
-                        // TODO: handle transformed mesh instances
-                default:
-                        Debug << "intersectable type not supported" << endl;
-                        break;
-                }
-                
-        if (mesh) // copy the mesh data to polygons
-                {
-                        mBox.Include(object->GetBox()); // add to BSP tree aabb
-                        AddMeshToPolygons(mesh, polys, mRootCell);
-                }
-        }
-
-        return (int)polys.size();
-}
-
-void VspBspTree::Construct(const VssRayContainer &sampleRays)
-{
-    mStat.nodes = 1;
-        mBox.Initialize();      // initialise BSP tree bounding box
-        
-        PolygonContainer polys;
-        RayInfoContainer *rays = new RayInfoContainer();
-
-        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
-
-        long startTime = GetTime();
-
-        Debug << "**** Extracting polygons from rays ****\n";
-
-        //-- extract polygons intersected by the rays
-        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
-        {
-                VssRay *ray = *rit;
-        
-                if (ray->mTerminationObject)
-                {
-                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
-                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
-                }
-
-                if (ray->mOriginObject)
-                {
-                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
-                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
-                }
-        }
-
-        //-- store rays
-        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
-        {
-                rays->push_back(RayInfo(*rit));
-        }
-
-        mStat.polys = (int)polys.size();
-
-        Debug << "**** Finished polygon extraction ****" << endl;
-        Debug << (int)polys.size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
-        Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
-
-        Construct(new PolygonContainer(polys), rays);
-
-        // clean up polygons because they are not deleted in the tree
-        CLEAR_CONTAINER(polys);
-}
-
-void VspBspTree::Construct(PolygonContainer *polys, RayInfoContainer *rays)
-{
-        std::stack<VspBspTraversalData> tStack;
-
-        mRoot = new VspBspLeaf();
-
-        VspBspNodeGeometry *cell = new VspBspNodeGeometry();
-        ConstructGeometry(mRoot, *cell);
-
-        VspBspTraversalData tData(mRoot, polys, 0, mRootCell, rays, 
-                                                   ComputePvsSize(*rays), cell->GetArea(), cell);
-
-        tStack.push(tData);
-
-        mStat.Start();
-        cout << "Contructing bsp tree ... ";
-        long startTime = GetTime();
-        while (!tStack.empty()) 
-        {
-                tData = tStack.top();
-
-            tStack.pop();
-
-                // subdivide leaf node
-                VspBspNode *r = Subdivide(tStack, tData);
-
-                if (r == mRoot)
-                        Debug << "BSP tree construction time spent at root: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
-        }
-
-        cout << "finished\n";
-
-        mStat.Stop();
-}
-
-bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
-{
-        return 
-                (((int)data.mRays->size() <= mTermMinRays) ||
-                 (data.mPvs <= mTermMinPvs) ||
-                 (data.mArea <= mTermMinArea) ||
-                 (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
-                 (data.mDepth >= mTermMaxDepth));
-}
-
-VspBspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, VspBspTraversalData &tData)
-{
-        //-- terminate traversal  
-        if (TerminationCriteriaMet(tData))              
-        {
-                VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(tData.mNode);
-        
-                BspViewCell *viewCell = new BspViewCell();
-                
-                leaf->SetViewCell(viewCell);
-                //TODO viewCell->mVspBspLeaves.push_back(leaf);
-
-                //-- add pvs
-                if (viewCell != mRootCell)
-                {
-                        int conSamp = 0, sampCon = 0;
-                        leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
-                        
-                        mStat.contributingSamples += conSamp;
-                        mStat.sampleContributions += sampCon;
-                }
-
-                EvaluateLeafStats(tData);
-                
-                //-- clean up
-
-                DEL_PTR(tData.mPolygons);
-                DEL_PTR(tData.mRays);
-                DEL_PTR(tData.mGeometry);
-                
-                return leaf;
-        }
-
-        //-- continue subdivision
-        PolygonContainer coincident;
-        
-        VspBspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
-                                                                new RayInfoContainer(), 0, 0, new VspBspNodeGeometry());
-        VspBspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell, 
-                                                           new RayInfoContainer(), 0, 0, new VspBspNodeGeometry());
-
-        // create new interior node and two leaf nodes
-        VspBspInterior *interior = 
-                SubdivideNode(tData, tFrontData, tBackData, coincident);
-
-#ifdef _DEBUG   
-        if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
-        {       for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
-                        Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
-                Debug << endl;}
-#endif
-
-        // don't need coincident polygons anymory
-        CLEAR_CONTAINER(coincident);
-
-        // push the children on the stack
-        tStack.push(tFrontData);
-        tStack.push(tBackData);
-
-        // cleanup
-        DEL_PTR(tData.mNode);
-        DEL_PTR(tData.mPolygons);
-        DEL_PTR(tData.mRays);
-        DEL_PTR(tData.mGeometry);               
-
-        return interior;
-}
-
-VspBspInterior *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
-                                                                        VspBspTraversalData &frontData,
-                                                                        VspBspTraversalData &backData,
-                                                                        PolygonContainer &coincident)
-{
-        mStat.nodes += 2;
-        
-        VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(tData.mNode);
-        // select subdivision plane
-        VspBspInterior *interior = 
-                new VspBspInterior(SelectPlane(leaf, tData)); 
-
-#ifdef _DEBUG
-        Debug << interior << endl;
-#endif
-        
-        // subdivide rays into front and back rays
-        SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);
-        
-        // subdivide polygons with plane
-        mStat.splits += interior->SplitPolygons(*tData.mPolygons, 
-                                                    *frontData.mPolygons, 
-                                                                                        *backData.mPolygons, 
-                                                                                        coincident);
-
-    // compute pvs
-        frontData.mPvs = ComputePvsSize(*frontData.mRays);
-        backData.mPvs = ComputePvsSize(*backData.mRays);
-
-        // split geometry and compute area
-        if (1)
-        {
-                tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
-                                                                           *backData.mGeometry, 
-                                                                           *this, 
-                                                                           interior->mPlane);
-        
-                
-                frontData.mArea = frontData.mGeometry->GetArea();
-                backData.mArea = backData.mGeometry->GetArea();
-        }
-
-        // compute accumulated ray length
-        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
-        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
-
-        //-- create front and back leaf
-
-        VspBspInterior *parent = leaf->GetParent();
-
-        // replace a link from node's parent
-        if (!leaf->IsRoot())
-        {
-                parent->ReplaceChildLink(leaf, interior);
-                interior->SetParent(parent);
-        }
-        else // new root
-        {
-                mRoot = interior;
-        }
-
-        // and setup child links
-        interior->SetupChildLinks(new VspBspLeaf(interior), new VspBspLeaf(interior));
-        
-        frontData.mNode = interior->mFront;
-        backData.mNode = interior->mBack;
-        
-        //DEL_PTR(leaf);
-        return interior;
-}
-
-void VspBspTree::SortSplitCandidates(const PolygonContainer &polys, 
-                                                                  const int axis, 
-                                                                  vector<SortableEntry> &splitCandidates) const
-{
-  splitCandidates.clear();
-  
-  int requestedSize = 2 * (int)polys.size();
-  // creates a sorted split candidates array  
-  splitCandidates.reserve(requestedSize);
-  
-  PolygonContainer::const_iterator it, it_end = polys.end();
-
-  AxisAlignedBox3 box;
-
-  // insert all queries 
-  for(it = polys.begin(); it != it_end; ++ it)
-  {
-          box.Initialize();
-          box.Include(*(*it));
-          
-          splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
-      splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
-  }
-  
-  stable_sort(splitCandidates.begin(), splitCandidates.end());
-}
-
-
-float VspBspTree::BestCostRatio(const PolygonContainer &polys,
-                                                         const AxisAlignedBox3 &box,
-                                                         const int axis,
-                                                         float &position,
-                                                         int &objectsBack,
-                                                         int &objectsFront) const
-{
-        vector<SortableEntry> splitCandidates;
-
-        SortSplitCandidates(polys, axis, splitCandidates);
-        
-        // go through the lists, count the number of objects left and right
-        // and evaluate the following cost funcion:
-        // C = ct_div_ci  + (ol + or)/queries
-        
-        int objectsLeft = 0, objectsRight = (int)polys.size();
-        
-        float minBox = box.Min(axis);
-        float maxBox = box.Max(axis);
-        float boxArea = box.SurfaceArea();
-  
-        float minBand = minBox + mSplitBorder * (maxBox - minBox);
-        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
-        
-        float minSum = 1e20f;
-        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
-
-        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci) 
-        {
-                switch ((*ci).type) 
-                {
-                        case SortableEntry::POLY_MIN:
-                                ++ objectsLeft;
-                                break;
-                        case SortableEntry::POLY_MAX:
-                            -- objectsRight;
-                                break;
-                        default:
-                                break;
-                }
-                
-                if ((*ci).value > minBand && (*ci).value < maxBand) 
-                {
-                        AxisAlignedBox3 lbox = box;
-                        AxisAlignedBox3 rbox = box;
-                        lbox.SetMax(axis, (*ci).value);
-                        rbox.SetMin(axis, (*ci).value);
-
-                        float sum = objectsLeft * lbox.SurfaceArea() + 
-                                                objectsRight * rbox.SurfaceArea();
-      
-                        if (sum < minSum) 
-                        {
-                                minSum = sum;
-                                position = (*ci).value;
-
-                                objectsBack = objectsLeft;
-                                objectsFront = objectsRight;
-                        }
-                }
-        }
-  
-        float oldCost = (float)polys.size();
-        float newCost = mAaCtDivCi + minSum / boxArea;
-        float ratio = newCost / oldCost;
-
-
-#if 0
-  Debug << "====================" << endl;
-  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
-      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
-#endif
-  return ratio;
-}
-
-bool VspBspTree::SelectAxisAlignedPlane(Plane3 &plane, 
-                                                                         const PolygonContainer &polys) const
-{
-        AxisAlignedBox3 box;
-        box.Initialize();
-        
-        // create bounding box of region
-        Polygon3::IncludeInBox(polys, box);
-        
-        int objectsBack = 0, objectsFront = 0;
-        int axis = 0;
-        float costRatio = MAX_FLOAT;
-        Vector3 position;
-
-        //-- area subdivision
-        for (int i = 0; i < 3; ++ i) 
-        {
-                float p = 0;
-                float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
-                
-                if (r < costRatio)
-                {
-                        costRatio = r;
-                        axis = i;
-                        position = p;
-                }
-        }
-        
-        if (costRatio >= mMaxCostRatio)
-                return false;
-
-        Vector3 norm(0,0,0); norm[axis] = 1.0f;
-        plane = Plane3(norm, position);
-
-        return true;
-}
-
-Plane3 VspBspTree::SelectPlane(VspBspLeaf *leaf, VspBspTraversalData &data)
-{
-        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
-                ((int)data.mRays->size() > mTermMinRaysForAxisAligned))
-        {
-                Plane3 plane;
-                if (SelectAxisAlignedPlane(plane, *data.mPolygons)) // TODO: should be rays
-                        return plane;
-        }
-
-        // simplest strategy: just take next polygon
-        if (mSplitPlaneStrategy & RANDOM_POLYGON)
-        {
-        if (!data.mPolygons->empty())
-                {
-                        Polygon3 *nextPoly = (*data.mPolygons)[Random((int)data.mPolygons->size())];
-                        return nextPoly->GetSupportingPlane();
-                }
-                else
-                {
-                        //-- choose plane on midpoint of a ray
-                        const int candidateIdx = Random((int)data.mRays->size());
-                                                                
-                        const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
-                        const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();
-
-                        const Vector3 pt = (maxPt + minPt) * 0.5;
-
-                        const Vector3 normal =  (*data.mRays)[candidateIdx].mRay->GetDir();
-                        
-                        return Plane3(normal, pt);
-                }
-
-                return Plane3();
-        }
-
-        // use heuristics to find appropriate plane
-        return SelectPlaneHeuristics(leaf, data); 
-}
-
-Plane3 VspBspTree::SelectPlaneHeuristics(VspBspLeaf *leaf, VspBspTraversalData &data)
-{
-        float lowestCost = MAX_FLOAT;
-        Plane3 bestPlane;
-        Plane3 plane;
-
-        int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
-        
-        int candidateIdx = limit;
-
-        for (int i = 0; i < limit; ++ i)
-        {
-                candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
-                
-                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
-                // evaluate current candidate
-                const float candidateCost = 
-                        SplitPlaneCost(poly->GetSupportingPlane(), data);
-
-                if (candidateCost < lowestCost)
-                {
-                        bestPlane = poly->GetSupportingPlane();
-                        lowestCost = candidateCost;
-                }
-        }
-        
-        //Debug << "lowest: " << lowestCost << endl;
-
-        //-- choose candidate planes extracted from rays
-        // we currently use different two methods chosen with
-        // equal probability
-        
-        RayInfoContainer *rays = data.mRays;
-        // take 3 ray endpoints, where two are minimum and one a maximum
-        // point or the other way round
-        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
-        {
-                candidateIdx = Random((int)rays->size());
-        
-                RayInfo rayInf = (*rays)[candidateIdx];
-                const Vector3 minPt = rayInf.ExtrapOrigin();
-                const Vector3 maxPt = rayInf.ExtrapTermination();
-
-                const Vector3 pt = (maxPt + minPt) * 0.5;
-                const Vector3 normal = Normalize(rayInf.mRay->GetDir());
-
-                plane = Plane3(normal, pt);
-        
-                const float candidateCost = SplitPlaneCost(plane, data);
-
-                if (candidateCost < lowestCost)
-                {
-                        bestPlane = plane;      
-                        lowestCost = candidateCost;
-                }
-        }
-
-        // take plane normal as plane normal and the midpoint of the ray.
-        // PROBLEM: does not resemble any point where visibility is likely to change
-        //Debug << "lowest: " << lowestCost << endl;
-        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
-        {
-                Vector3 pt[3];
-                int idx[3];
-                int cmaxT = 0;
-                int cminT = 0;
-                bool chooseMin = false;
-
-                for (int j = 0; j < 3; j ++)
-                {
-                        idx[j] = Random((int)rays->size() * 2);
-                                
-                        if (idx[j] >= (int)rays->size())
-                        {
-                                idx[j] -= (int)rays->size();
-                                
-                                chooseMin = (cminT < 2);
-                        }
-                        else
-                                chooseMin = (cmaxT < 2);
-
-                        RayInfo rayInf = (*rays)[idx[j]];
-                        pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
-                }       
-                        
-                plane = Plane3(pt[0], pt[1], pt[2]);
-
-                const float candidateCost = SplitPlaneCost(plane, data);
-
-                if (candidateCost < lowestCost)
-                {
-                        //Debug << "choose ray plane 2: " << candidateCost << endl;
-                        bestPlane = plane;
-                        
-                        lowestCost = candidateCost;
-                }
-        }       
-
-#ifdef _DEBUG
-        Debug << "plane lowest cost: " << lowestCost << endl;
-#endif
-        return bestPlane;
-}
-
-int VspBspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
-{
-        const int candidateIdx = Random(currentIdx --);
-
-        // swap candidates to avoid testing same plane 2 times
-        std::swap(polys[currentIdx], polys[candidateIdx]);
-        
-        return currentIdx;
-        //return Random((int)polys.size());
-}
-
-
-bool VspBspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
-{
-        maxT = 1e6;
-        minT = 0;
-        
-        // test with tree bounding box
-        if (!mBox.GetMinMaxT(ray, &minT, &maxT))
-                return false;
-
-        if (minT < 0) // start ray from origin
-                minT = 0;
-
-        // bound ray or line segment
-        if ((ray.GetType() == Ray::LOCAL_RAY) && 
-            !ray.intersections.empty() && 
-                (ray.intersections[0].mT <= maxT))
-        {
-                maxT = ray.intersections[0].mT;
-        }
-
-        return true;
-}
-
-inline void VspBspTree::GenerateUniqueIdsForPvs()
-{
-        Intersectable::NewMail(); sBackId = ViewCell::sMailId;
-        Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
-        Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
-}
-
-float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
-                                                                 const VspBspTraversalData &data)
-{
-        float val = 0;
-
-        float sumBalancedRays = 0;
-        float sumRaySplits = 0;
-        
-        int backId = 0;
-        int frontId = 0;
-        int frontAndBackId = 0;
-
-        int frontPvs = 0;
-        int backPvs = 0;
-
-        // probability that view point lies in child
-        float pOverall = 0;
-        float pFront = 0;
-        float pBack = 0;
-
-        const bool pvsUseLen = false;
-
-        if (mSplitPlaneStrategy & PVS)
-        {
-                // create unique ids for pvs heuristics
-                GenerateUniqueIdsForPvs();
-        
-                if (mPvsUseArea) // use front and back cell areas to approximate volume
-                {       
-                        // construct child geometry with regard to the candidate split plane
-                        VspBspNodeGeometry frontCell;
-                        VspBspNodeGeometry backCell;
-
-                        data.mGeometry->SplitGeometry(frontCell, backCell, *this, candidatePlane);
-                
-                        pFront = frontCell.GetArea();
-                        pBack = backCell.GetArea();
-
-                        pOverall = data.mArea;
-                }
-        }
-                        
-        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
-
-        for (rit = data.mRays->begin(); rit != data.mRays->end(); ++ rit)
-        {
-                VssRay *ray = (*rit).mRay;
-                const int cf = (*rit).ComputeRayIntersection(candidatePlane, (*rit).mRay->mT);
-
-                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-                {
-                        sumBalancedRays += cf;
-                }
-                
-                if (mSplitPlaneStrategy & BALANCED_RAYS)
-                {
-                        if (cf == 0)
-                                ++ sumRaySplits;
-                }
-
-                if (mSplitPlaneStrategy & PVS)
-                {
-                        // in case the ray intersects an object
-                        // assure that we only count the object 
-                        // once for the front and once for the back side of the plane
-                        
-                        // add the termination object
-                        AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
-                        
-                        // add the source object
-                        AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
-                        
-                        // use number or length of rays to approximate volume
-                        if (!mPvsUseArea)
-                        {
-                                float len = 1;
-
-                                if (pvsUseLen) // use length of rays
-                                        len = (*rit).SqrSegmentLength();
-                        
-                                pOverall += len;
-
-                                if (cf == 1)
-                                        pFront += len;
-                                if (cf == -1)
-                                        pBack += len;
-                                if (cf == 0)
-                                {
-                                        // use length of rays to approximate volume
-                                        if (pvsUseLen) 
-                                        {
-                                                float newLen = len * 
-                                                        ((*rit).GetMaxT() - (*rit).mRay->mT) / 
-                                                        ((*rit).GetMaxT() - (*rit).GetMinT());
-                
-                                                if (candidatePlane.Side((*rit).ExtrapOrigin()) <= 0)
-                                                {
-                                                        pBack += newLen;
-                                                        pFront += len - newLen;
-                                                }
-                                                else
-                                                {
-                                                        pFront += newLen;
-                                                        pBack += len - newLen;
-                                                }
-                                        }
-                                        else
-                                        {
-                                                ++ pFront;
-                                                ++ pBack;
-                                        }
-                                }
-                        }
-                }
-        }
-
-        const float raysSize = (float)data.mRays->size() + Limits::Small;
-
-        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
-                val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
-
-        if (mSplitPlaneStrategy & BALANCED_RAYS)
-                val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
-
-        const float denom = pOverall * (float)data.mPvs * 2.0f + Limits::Small;
-
-        if (mSplitPlaneStrategy & PVS)
-        {
-                val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
-
-                // give penalty to unbalanced split
-                if (0)
-                if (((pFront * 0.2 + Limits::Small) > pBack) || 
-                        (pFront < (pBack * 0.2 + Limits::Small)))
-                        val += 0.5;
-        }
-
-#ifdef _DEBUG
-        Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
-                  << " frontpvs: " << frontPvs << " pFront: " << pFront 
-                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
-#endif
-        return val;
-}
-
-void VspBspTree::AddObjToPvs(Intersectable *obj,
-                                                  const int cf,
-                                                  int &frontPvs,
-                                                  int &backPvs) const
-{
-        if (!obj)
-                return;
-        // TODO: does this really belong to no pvs?
-        //if (cf == Ray::COINCIDENT) return;
-
-        // object belongs to both PVS
-        if (cf >= 0)
-        {
-                if ((obj->mMailbox != sFrontId) && 
-                        (obj->mMailbox != sFrontAndBackId))
-                {
-                        ++ frontPvs;
-
-                        if (obj->mMailbox == sBackId)
-                                obj->mMailbox = sFrontAndBackId;        
-                        else
-                                obj->mMailbox = sFrontId;                                                               
-                }
-        }
-        
-        if (cf <= 0)
-        {
-                if ((obj->mMailbox != sBackId) &&
-                        (obj->mMailbox != sFrontAndBackId))
-                {
-                        ++ backPvs;
-
-                        if (obj->mMailbox == sFrontId)
-                                obj->mMailbox = sFrontAndBackId; 
-                        else
-                                obj->mMailbox = sBackId;                                
-                }
-        }
-}
-
-void VspBspTree::CollectLeaves(vector<VspBspLeaf *> &leaves) const
-{
-        stack<VspBspNode *> nodeStack;
-        nodeStack.push(mRoot);
- 
-        while (!nodeStack.empty()) 
-        {
-                VspBspNode *node = nodeStack.top();
-    
-                nodeStack.pop();
-    
-                if (node->IsLeaf()) 
-                {
-                        VspBspLeaf *leaf = (VspBspLeaf *)node;          
-                        leaves.push_back(leaf);
-                } 
-                else 
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-                        nodeStack.push(interior->GetBack());
-                        nodeStack.push(interior->GetFront());
-                }
-        }
-}
-
-AxisAlignedBox3 VspBspTree::GetBoundingBox() const
-{
-        return mBox;
-}
-
-VspBspNode *VspBspTree::GetRoot() const
-{
-        return mRoot;
-}
-
-void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
-{
-        // the node became a leaf -> evaluate stats for leafs
-        VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(data.mNode);
-
-        // store maximal and minimal depth
-        if (data.mDepth > mStat.maxDepth)
-                mStat.maxDepth = data.mDepth;
-
-        if (data.mDepth < mStat.minDepth)
-                mStat.minDepth = data.mDepth;
-
-        // accumulate depth to compute average depth
-        mStat.accumDepth += data.mDepth;
-        
-
-        if (data.mDepth >= mTermMaxDepth)
-                ++ mStat.maxDepthNodes;
-
-        if (data.mPvs < mTermMinPvs)
-                ++ mStat.minPvsNodes;
-
-        if ((int)data.mRays->size() < mTermMinRays)
-                ++ mStat.minRaysNodes;
-
-        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
-                ++ mStat.maxRayContribNodes;
-        
-        if (data.mGeometry->GetArea() <= mTermMinArea) 
-                ++ mStat.minAreaNodes;
-
-#ifdef _DEBUG
-        Debug << "BSP stats: "
-                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
-                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
-                  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
-                  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << mTermMinPolys << "), "
-                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), " 
-                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
-                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
-#endif
-}
-
-int VspBspTree::CastRay(Ray &ray)
-{
-        int hits = 0;
-  
-        stack<VspBspRayTraversalData> tStack;
-  
-        float maxt, mint;
-
-        if (!BoundRay(ray, mint, maxt))
-                return 0;
-
-        Intersectable::NewMail();
-
-        Vector3 entp = ray.Extrap(mint);
-        Vector3 extp = ray.Extrap(maxt);
-  
-        VspBspNode *node = mRoot;
-        VspBspNode *farChild = NULL;
-        
-        while (1)
-        {
-                if (!node->IsLeaf()) 
-                {
-                        VspBspInterior *in = (VspBspInterior *) node;
-                        
-                        Plane3 *splitPlane = in->GetPlane();
-
-                        int entSide = splitPlane->Side(entp);
-                        int extSide = splitPlane->Side(extp);
-
-                        Vector3 intersection;
-
-                        if (entSide < 0)
-                        {
-                                node = in->GetBack();
-
-                                if(extSide <= 0) // plane does not split ray => no far child
-                                        continue;
-                                        
-                                farChild = in->GetFront(); // plane splits ray
-
-                        } else if (entSide > 0)
-                        {
-                                node = in->GetFront();
-
-                                if (extSide >= 0) // plane does not split ray => no far child
-                                        continue;
-
-                                farChild = in->GetBack(); // plane splits ray                   
-                        }
-                        else // ray and plane are coincident
-                        {
-                                // WHAT TO DO IN THIS CASE ?
-                                //break;
-                                node = in->GetFront();
-                                continue;
-                        }
-
-                        // push data for far child
-                        tStack.push(VspBspRayTraversalData(farChild, extp, maxt));
-
-                        // find intersection of ray segment with plane
-                        float t;
-                        extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &t);
-                        maxt *= t;
-                        
-                } else // reached leaf => intersection with view cell
-                {
-                        VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(node);
-      
-                        if (!leaf->mViewCell->Mailed())
-                        {
-                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
-                                leaf->mViewCell->Mail();
-                                ++ hits;
-                        }
-                        
-                        //-- fetch the next far child from the stack
-                        if (tStack.empty())
-                                break;
-      
-                        entp = extp;
-                        mint = maxt; // NOTE: need this?
-
-                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
-                                break;
-
-                        VspBspRayTraversalData &s = tStack.top();
-
-                        node = s.mNode;
-                        extp = s.mExitPoint;
-                        maxt = s.mMaxT;
-
-                        tStack.pop();
-                }
-        }
-
-        return hits;
-}
-
-bool VspBspTree::Export(const string filename)
-{
-        Exporter *exporter = Exporter::GetExporter(filename);
-
-        if (exporter) 
-        {
-                //exporter->ExportVspBspTree(*this);
-                return true;
-        }       
-
-        return false;
-}
-
-void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
-{
-        stack<VspBspNode *> nodeStack;
-        nodeStack.push(mRoot);
-
-        ViewCell::NewMail();
-
-        while (!nodeStack.empty()) 
-        {
-                VspBspNode *node = nodeStack.top();
-                nodeStack.pop();
-
-                if (node->IsLeaf()) 
-                {
-                        ViewCell *viewCell = dynamic_cast<VspBspLeaf *>(node)->mViewCell;
-
-                        if (!viewCell->Mailed()) 
-                        {
-                                viewCell->Mail();
-                                viewCells.push_back(viewCell);
-                        }
-                }
-                else 
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-                        nodeStack.push(interior->mFront);
-                        nodeStack.push(interior->mBack);
-                }
-        }
-}
-
-void VspBspTree::EvaluateViewCellsStats(VspBspViewCellsStatistics &stat) const
-{
-        stat.Reset();
-
-        stack<VspBspNode *> nodeStack;
-        nodeStack.push(mRoot);
-
-        ViewCell::NewMail();
-
-        // exclude root cell
-        mRootCell->Mail();
-
-        while (!nodeStack.empty()) 
-        {
-                VspBspNode *node = nodeStack.top();
-                nodeStack.pop();
-
-                if (node->IsLeaf()) 
-                {
-                        ++ stat.bspLeaves;
-
-                        BspViewCell *viewCell = dynamic_cast<VspBspLeaf *>(node)->mViewCell;
-
-                        if (!viewCell->Mailed()) 
-                        {
-                                viewCell->Mail();
-                                
-                                ++ stat.viewCells;
-                                const int pvsSize = viewCell->GetPvs().GetSize();
-
-                stat.pvs += pvsSize;
-
-                                if (pvsSize < 1)
-                                        ++ stat.emptyPvs;
-
-                                if (pvsSize > stat.maxPvs)
-                                        stat.maxPvs = pvsSize;
-
-                                if (pvsSize < stat.minPvs)
-                                        stat.minPvs = pvsSize;
-
-                                //TODO
-                                //if ((int)viewCell->mVspBspLeaves.size() > stat.maxVspBspLeaves)
-                                //      stat.maxVspBspLeaves = (int)viewCell->mVspBspLeaves.size();             
-                        }
-                }
-                else 
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-                        nodeStack.push(interior->mFront);
-                        nodeStack.push(interior->mBack);
-                }
-        }
-}
-
-VspBspTreeStatistics &VspBspTree::GetStat()
-{
-        return mStat;
-}
-
-float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
-{
-        float len = 0;
-
-        RayInfoContainer::const_iterator it, it_end = rays.end();
-
-        for (it = rays.begin(); it != it_end; ++ it)
-                len += (*it).SegmentLength();
-
-        return len;
-}
-
-int VspBspTree::SplitRays(const Plane3 &plane,
-                                           RayInfoContainer &rays, 
-                                           RayInfoContainer &frontRays, 
-                                           RayInfoContainer &backRays)
-{
-        int splits = 0;
-
-        while (!rays.empty())
-        {
-                RayInfo bRay = rays.back();
-                VssRay *ray = bRay.mRay;
-
-                rays.pop_back();
-        
-                // get classification and receive new t
-                const int cf =  bRay.ComputeRayIntersection(plane, ray->mT);
-
-                switch(cf)
-                {
-                case -1:
-                        backRays.push_back(bRay);
-                        break;
-                case 1:
-                        break;
-                        frontRays.push_back(bRay);
-                case 0: 
-                        //-- split ray
-
-                        // if start point behind plane
-                        if (plane.Side(bRay.ExtrapOrigin()) <= 0)
-                        {
-                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
-                                frontRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
-                        }
-                        else
-                        {
-                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
-                                backRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
-                        }
-                        break;
-                default:
-                        Debug << "Should not come here 4" << endl;
-                        break;
-                }
-        }
-
-        return splits;
-}
-
-void VspBspTree::ExtractHalfSpaces(VspBspNode *n, vector<Plane3> &halfSpaces) const
-{
-        VspBspNode *lastNode;
-        do
-        {
-                lastNode = n;
-
-                // want to get planes defining geometry of this node => don't take
-                // split plane of node itself
-                n = n->GetParent();
-                
-                if (n)
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(n);
-                        Plane3 halfSpace = *dynamic_cast<VspBspInterior *>(interior)->GetPlane();
-
-            if (interior->mFront != lastNode)
-                                halfSpace.ReverseOrientation();
-
-                        halfSpaces.push_back(halfSpace);
-                }
-        }
-        while (n);
-}
-
-void VspBspTree::ConstructGeometry(VspBspNode *n, VspBspNodeGeometry &cell) const
-{
-        PolygonContainer polys;
-        ConstructGeometry(n, polys);
-        cell.mPolys = polys;
-}
-
-void VspBspTree::ConstructGeometry(VspBspNode *n, PolygonContainer &cell) const
-{
-        vector<Plane3> halfSpaces;
-        ExtractHalfSpaces(n, halfSpaces);
-
-        PolygonContainer candidatePolys;
-
-        // bounded planes are added to the polygons (reverse polygons
-        // as they have to be outfacing
-        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
-        {
-                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
-                
-                if (p->Valid())
-                {
-                        candidatePolys.push_back(p->CreateReversePolygon());
-                        DEL_PTR(p);
-                }
-        }
-
-        // add faces of bounding box (also could be faces of the cell)
-        for (int i = 0; i < 6; ++ i)
-        {
-                VertexContainer vertices;
-        
-                for (int j = 0; j < 4; ++ j)
-                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
-
-                candidatePolys.push_back(new Polygon3(vertices));
-        }
-
-        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
-        {
-                // polygon is split by all other planes
-                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
-                {
-                        if (i == j) // polygon and plane are coincident
-                                continue;
-
-                        VertexContainer splitPts;
-                        Polygon3 *frontPoly, *backPoly;
-
-                        const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
-                        
-                        switch (cf)
-                        {
-                                case Polygon3::SPLIT:
-                                        frontPoly = new Polygon3();
-                                        backPoly = new Polygon3();
-
-                                        candidatePolys[i]->Split(halfSpaces[j], 
-                                                                                         *frontPoly, 
-                                                                                         *backPoly);
-
-                                        DEL_PTR(candidatePolys[i]);
-
-                                        if (frontPoly->Valid())
-                                                candidatePolys[i] = frontPoly;
-                                        else
-                                                DEL_PTR(frontPoly);
-
-                                        DEL_PTR(backPoly);
-                                        break;
-                                case Polygon3::BACK_SIDE:
-                                        DEL_PTR(candidatePolys[i]);
-                                        break;
-                                // just take polygon as it is
-                                case Polygon3::FRONT_SIDE:
-                                case Polygon3::COINCIDENT:
-                                default:
-                                        break;
-                        }
-                }
-                
-                if (candidatePolys[i])
-                        cell.push_back(candidatePolys[i]);
-        }
-}
-
-void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &cell) const
-{
-        /*
-        TODO
-        vector<VspBspLeaf *> leaves = vc->mBspLeaves;
-
-        vector<VspBspLeaf *>::const_iterator it, it_end = leaves.end();
-
-        for (it = leaves.begin(); it != it_end; ++ it)
-                ConstructGeometry(*it, vcGeom);*/
-}
-
-int VspBspTree::FindNeighbors(VspBspNode *n, vector<VspBspLeaf *> &neighbors, 
-                                                   const bool onlyUnmailed) const
-{
-        PolygonContainer cell;
-
-        ConstructGeometry(n, cell);
-
-        stack<VspBspNode *> nodeStack;
-        nodeStack.push(mRoot);
-                
-        // planes needed to verify that we found neighbor leaf.
-        vector<Plane3> halfSpaces;
-        ExtractHalfSpaces(n, halfSpaces);
-
-        while (!nodeStack.empty()) 
-        {
-                VspBspNode *node = nodeStack.top();
-                nodeStack.pop();
-
-                if (node->IsLeaf())
-                {
-            if (node != n && (!onlyUnmailed || !node->Mailed()))
-                        {
-                                // test all planes of current node if candidate really
-                                // is neighbour
-                                PolygonContainer neighborCandidate;
-                                ConstructGeometry(node, neighborCandidate);
-                                
-                                bool isAdjacent = true;
-                                for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
-                                {
-                                        const int cf = 
-                                                Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);
-
-                                        if (cf == Polygon3::BACK_SIDE)
-                                                isAdjacent = false;
-                                }
-
-                                if (isAdjacent)
-                                        neighbors.push_back(dynamic_cast<VspBspLeaf *>(node));
-
-                                CLEAR_CONTAINER(neighborCandidate);
-                        }
-                }
-                else 
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-        
-                        const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);
-
-                        if (cf == Polygon3::FRONT_SIDE)
-                                nodeStack.push(interior->mFront);
-                        else
-                                if (cf == Polygon3::BACK_SIDE)
-                                        nodeStack.push(interior->mBack);
-                                else 
-                                {
-                                        // random decision
-                                        nodeStack.push(interior->mBack);
-                                        nodeStack.push(interior->mFront);
-                                }
-                }
-        }
-        
-        CLEAR_CONTAINER(cell);
-        return (int)neighbors.size();
-}
-
-VspBspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
-{
-    stack<VspBspNode *> nodeStack;
-        nodeStack.push(mRoot);
-        
-        int mask = rand();
-  
-        while (!nodeStack.empty()) 
-        {
-                VspBspNode *node = nodeStack.top();
-                nodeStack.pop();
-          
-                if (node->IsLeaf()) 
-                {
-                        return dynamic_cast<VspBspLeaf *>(node);
-                } 
-                else 
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-                        
-                        VspBspNode *next;
-        
-                        PolygonContainer cell;
-
-                        // todo: not very efficient: constructs full cell everytime
-                        ConstructGeometry(interior, cell);
-
-                        const int cf = Polygon3::ClassifyPlane(cell, halfspace);
-
-                        if (cf == Polygon3::BACK_SIDE)
-                                next = interior->mFront;
-                        else
-                                if (cf == Polygon3::FRONT_SIDE)
-                                        next = interior->mFront;
-                        else 
-                        {
-                                // random decision
-                                if (mask & 1)
-                                        next = interior->mBack;
-                                else
-                                        next = interior->mFront;
-                                mask = mask >> 1;
-                        }
-
-                        nodeStack.push(next);
-                }
-        }
-        
-        return NULL;
-}
-
-VspBspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
-{
-        stack<VspBspNode *> nodeStack;
-        
-        nodeStack.push(mRoot);
-
-        int mask = rand();
-        
-        while (!nodeStack.empty()) 
-        {
-                VspBspNode *node = nodeStack.top();
-                nodeStack.pop();
-                
-                if (node->IsLeaf()) 
-                {
-                        if ( (!onlyUnmailed || !node->Mailed()) )
-                                return dynamic_cast<VspBspLeaf *>(node);
-                }
-                else 
-                {
-                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
-
-                        // random decision
-                        if (mask & 1)
-                                nodeStack.push(interior->mBack);
-                        else
-                                nodeStack.push(interior->mFront);
-
-                        mask = mask >> 1;
-                }
-        }
-        
-        return NULL;
-}
-
-int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
-{
-        int pvsSize = 0;
-
-        RayInfoContainer::const_iterator rit, rit_end = rays.end();
-
-        Intersectable::NewMail();
-
-        for (rit = rays.begin(); rit != rays.end(); ++ rit)
-        {
-                VssRay *ray = (*rit).mRay;
-                
-                if (ray->mOriginObject)
-                {
-                        if (!ray->mOriginObject->Mailed())
-                        {
-                                ray->mOriginObject->Mail();
-                                ++ pvsSize;
-                        }
-                }
-                if (ray->mTerminationObject)
-                {
-                        if (!ray->mTerminationObject->Mailed())
-                        {
-                                ray->mTerminationObject->Mail();
-                                ++ pvsSize;
-                        }
-                }
-        }
-
-        return pvsSize;
-}
-
-/*************************************************************
- *            VspBspNodeGeometry Implementation                 *
- *************************************************************/
-
-VspBspNodeGeometry::~VspBspNodeGeometry()
-{
-        CLEAR_CONTAINER(mPolys);
-}
-
-float VspBspNodeGeometry::GetArea() const 
-{ 
-        return Polygon3::GetArea(mPolys);
-}
-
-void VspBspNodeGeometry::SplitGeometry(VspBspNodeGeometry &front,
-                                                                        VspBspNodeGeometry &back,
-                                                                        const VspBspTree &tree,                                          
-                                                                        const Plane3 &splitPlane) const
-{       
-        // get cross section of new polygon
-        Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);
-
-        planePoly = SplitPolygon(planePoly, tree);
-
-        //-- plane poly splits all other cell polygons
-        for (int i = 0; i < (int)mPolys.size(); ++ i)
-        {
-                const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
-                        
-                // split new polygon with all previous planes
-                switch (cf)
-                {
-                        case Polygon3::SPLIT:
-                                {
-                                        Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);
-
-                                        Polygon3 *frontPoly = new Polygon3();
-                                        Polygon3 *backPoly = new Polygon3();
-                                
-                                        poly->Split(splitPlane, *frontPoly, *backPoly);
-
-                                        DEL_PTR(poly);
-
-                                        if (frontPoly->Valid())
-                                                front.mPolys.push_back(frontPoly);
-                                        if (backPoly->Valid())
-                                                back.mPolys.push_back(backPoly);
-                                }
-                                
-                                break;
-                        case Polygon3::BACK_SIDE:
-                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));                      
-                                break;
-                        case Polygon3::FRONT_SIDE:
-                                front.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));     
-                                break;
-                        case Polygon3::COINCIDENT:
-                                //front.mPolys.push_back(CreateReversePolygon(mPolys[i]));
-                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));
-                                break;
-                        default:
-                                break;
-                }
-        }
-
-        //-- finally add the new polygon to the child cells
-        if (planePoly)
-        {
-                // add polygon with normal pointing into positive half space to back cell
-                back.mPolys.push_back(planePoly);
-                // add polygon with reverse orientation to front cell
-                front.mPolys.push_back(planePoly->CreateReversePolygon());
-        }
-
-        //Debug << "returning new geometry " << mPolys.size() << " f: " << front.mPolys.size() << " b: " << back.mPolys.size() << endl;
-        //Debug << "old area " << GetArea() << " f: " << front.GetArea() << " b: " << back.GetArea() << endl;
-}
-
-Polygon3 *VspBspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
-                                                                                const VspBspTree &tree) const
-{
-        // polygon is split by all other planes
-        for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
-        {
-                Plane3 plane = mPolys[i]->GetSupportingPlane();
-
-                const int cf = 
-                        planePoly->ClassifyPlane(plane, 0.00001f);
-                        
-                // split new polygon with all previous planes
-                switch (cf)
-                {
-                        case Polygon3::SPLIT:
-                                {
-                                        Polygon3 *frontPoly = new Polygon3();
-                                        Polygon3 *backPoly = new Polygon3();
-
-                                        planePoly->Split(plane, *frontPoly, *backPoly);
-                                        DEL_PTR(planePoly);
-
-                                        if (backPoly->Valid())
-                                                planePoly = backPoly;
-                                        else
-                                                DEL_PTR(backPoly);
-                                }
-                                break;
-                        case Polygon3::FRONT_SIDE:
-                                DEL_PTR(planePoly);
-                break;
-                        // polygon is taken as it is
-                        case Polygon3::BACK_SIDE:
-                        case Polygon3::COINCIDENT:
-                        default:
-                                break;
-                }
-        }
-
-        return planePoly;
-}
-
-void VspBspViewCellsStatistics::Print(ostream &app) const
-{
-        app << "===== VspBspViewCells statistics ===============\n";
-
-        app << setprecision(4);
-
-        //app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
-
-        app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
-
-        app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
-
-        app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
-
-        app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
-
-        app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
-
-        app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
-
-        app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgVspBspLeaves() << endl;
-
-        app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxVspBspLeaves << endl;
-        
-        app << "===== END OF VspBspViewCells statistics ==========\n";
-}
-
-BspViewCell *VspBspTree::GetRootCell() const
-{
-        return mRootCell;
-}

trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.h

@@ -9 +9 @@
 #include "VssRay.h"
 #include "RayInfo.h"
+#include "ViewCellBsp.h"
 
 class ViewCell;
@@ -14 +15 @@
 class Plane3;
 class VspBspTree;  
-class VspBspInterior;
-class VspBspNode;
+class BspInterior;
+class BspNode;
 class AxisAlignedBox3;
 class Ray;
 
-class VspBspNodeGeometry
-{
-public:
-        VspBspNodeGeometry()
-        {};  
-
-        ~VspBspNodeGeometry();
-
-        float GetArea() const;
-        
-        /** Computes new cell based on the old cell definition and a new split plane
-                @param side indicates which side of the halfspace 
-        */
-        void SplitGeometry(VspBspNodeGeometry &front, 
-                                           VspBspNodeGeometry &back,
-                                           const VspBspTree &tree,
-                                           const Plane3 &splitPlane) const;
-
-        Polygon3 *SplitPolygon(Polygon3 *poly, const VspBspTree &tree) const;
-
-        PolygonContainer mPolys;
-};
-
-/** Data structure used for optimized ray casting.
+/*class BspNodeGeometry;
+class BspTreeStatistics;
+class BspViewCellsStatistics;
+class BspNode;
+class BspLeaf;
+class BspInterior;
 */
-struct VspBspRayTraversalData 
-{
-    VspBspNode *mNode;
-    Vector3 mExitPoint;
-    float mMaxT;
-    
-    VspBspRayTraversalData() {}
-
-    VspBspRayTraversalData(VspBspNode *n, const Vector3 &extp, const float maxt):
-    mNode(n), mExitPoint(extp), mMaxT(maxt) 
-        {}
-};
-
-class VspBspTreeStatistics: public StatisticsBase
-{
-public:
-        // total number of nodes
-        int nodes;
-        // number of splits
-        int splits;
-        // totals number of rays
-        int rays;
-        // maximal reached depth
-        int maxDepth;
-        // minimal depth
-        int minDepth;
-        
-        // max depth nodes
-        int maxDepthNodes;
-        // minimum depth nodes
-        int minDepthNodes;
-        // max depth nodes
-        int minPvsNodes;
-        // nodes with minimum PVS
-        int minRaysNodes;
-        // max ray contribution nodes
-        int maxRayContribNodes;
-        // minimum area nodes
-        int minAreaNodes;
-
-        // max number of rays per node
-        int maxObjectRefs;
-        // accumulated depth (used to compute average)
-        int accumDepth;
-        // number of initial polygons
-        int polys;
-        /// samples contributing to pvs
-        int contributingSamples;
-        /// sample contributions to pvs
-        int sampleContributions;
-        /// largest pvs
-        int largestPvs;
-
-        // Constructor
-        VspBspTreeStatistics() 
-        {
-                Reset();
-        }
-
-        int Nodes() const {return nodes;}
-        int Interior() const { return nodes / 2; }
-        int Leaves() const { return (nodes / 2) + 1; }
-        
-        // TODO: computation wrong
-        double AvgDepth() const { return accumDepth / (double)Leaves();}; 
-  
-        void Reset() 
-        {
-                nodes = 0;
-                splits = 0;
-                
-                maxDepth = 0;
-                minDepth = 99999;
-                polys = 0;
-                accumDepth = 0;
-        
-                maxDepthNodes = 0;
-                minPvsNodes = 0;
-                minRaysNodes = 0;
-                maxRayContribNodes = 0;
-                minAreaNodes = 0;
-
-                contributingSamples = 0;
-                sampleContributions = 0;
-        }
-
-        void Print(ostream &app) const;
-
-        friend ostream &operator<<(ostream &s, const VspBspTreeStatistics &stat) 
-        {
-                stat.Print(s);
-                return s;
-        } 
-};
-
-class VspBspViewCellsStatistics: public StatisticsBase
-{
-public:
-
-        /// number of view cells
-        int viewCells;
-
-        /// size of the PVS
-        int pvs;
-  
-        /// largest PVS of all view cells
-        int maxPvs;
-
-        /// smallest PVS of all view cells
-        int minPvs;
-
-        /// view cells with empty PVS
-        int emptyPvs;
-
-        /// number of bsp leaves covering the view space
-        int bspLeaves;
-
-        /// largest number of leaves covered by one view cell  
-        int maxVspBspLeaves;
-
-    // Constructor
-        VspBspViewCellsStatistics() 
-        {
-                Reset();
-        }
-
-        double AvgVspBspLeaves() const {return (double)bspLeaves / (double)viewCells;};
-        double AvgPvs() const {return (double)pvs / (double)viewCells;};
-  
-        void Reset() 
-        {
-                viewCells = 0;  
-                pvs = 0;  
-                maxPvs = 0;
-
-                minPvs = 999999;
-                emptyPvs = 0;
-                bspLeaves = 0;
-                maxVspBspLeaves = 0;
-        }
-
-        void Print(ostream &app) const;
-
-        friend ostream &operator<<(ostream &s, const VspBspViewCellsStatistics &stat) 
-        {
-                stat.Print(s);
-                return s;
-        }  
-};
-
-/**
-    VspBspNode abstract class serving for interior and leaf node implementation
-*/
-class VspBspNode 
-{
-        friend class VspBspTree;
-
-public:
-        VspBspNode();
-        virtual ~VspBspNode(){};
-        VspBspNode(VspBspInterior *parent);
-
-        /** Determines whether this node is a leaf or not
-        @return true if leaf
-        */
-        virtual bool IsLeaf() const = 0;
-
-        /** Determines whether this node is a root
-        @return true if root
-        */
-        virtual bool IsRoot() const;
-
-        /** Returns parent node.
-        */
-        VspBspInterior *GetParent();
-
-        /** Sets parent node.
-        */
-        void SetParent(VspBspInterior *parent);
-
-        
-        static int sMailId;
-        int mMailbox;
-  
-        void Mail() { mMailbox = sMailId; }
-        static void NewMail() { ++ sMailId; }
-        bool Mailed() const { return mMailbox == sMailId; }
-
-protected:
-
-        /// parent of this node
-        VspBspInterior *mParent;
-};
-
-/** BSP interior node implementation 
-*/
-class VspBspInterior : public VspBspNode 
-{
-        friend class VspBspTree;
-public:
-        /** Standard contructor taking split plane as argument.
-        */
-        VspBspInterior(const Plane3 &plane);
-        ~VspBspInterior();
-        /** @return false since it is an interior node 
-        */
-        bool IsLeaf() const;
-
-        VspBspNode *GetBack();
-        VspBspNode *GetFront();
-
-        Plane3 *GetPlane();
-
-        void ReplaceChildLink(VspBspNode *oldChild, VspBspNode *newChild);
-        void SetupChildLinks(VspBspNode *b, VspBspNode *f);
-
-        /** Splits polygons with respect to the split plane.
-                @param polys the polygons to be split. the polygons are consumed and
-                           distributed to the containers frontPolys, backPolys, coincident.
-                @param frontPolys returns the polygons in the front of the split plane
-                @param backPolys returns the polygons in the back of the split plane
-                @param coincident returns the polygons coincident to the split plane
-
-                @returns the number of splits   
-        */
-        int SplitPolygons(PolygonContainer &polys, 
-                                          PolygonContainer &frontPolys, 
-                                          PolygonContainer &backPolys, 
-                                          PolygonContainer &coincident);
-
-        friend ostream &operator<<(ostream &s, const VspBspInterior &A)
-        {
-                return s << A.mPlane;
-        }
-
-protected:
-
-        /// Splitting plane corresponding to this node
-        Plane3 mPlane;
-        /// back node
-        VspBspNode *mBack;
-        /// front node
-        VspBspNode *mFront;
-};
-
-/** BSP leaf node implementation.
-*/
-class VspBspLeaf : public VspBspNode 
-{
-        friend class VspBspTree;
-
-public:
-        VspBspLeaf();
-        VspBspLeaf(BspViewCell *viewCell);
-        VspBspLeaf(VspBspInterior *parent);
-        VspBspLeaf(VspBspInterior *parent, BspViewCell *viewCell);
-
-        /** @return true since it is an interior node 
-        */
-        bool IsLeaf() const;
-        
-        /** Returns pointer of view cell.
-        */
-        BspViewCell *GetViewCell() const;
-
-        /** Sets pointer to view cell.
-        */
-        void SetViewCell(BspViewCell *viewCell);
-
-        /** Adds ray sample contributions to the PVS.
-                @param sampleContributions the number contributions of the samples
-                @param contributingSampels the number of contributing rays
-                
-        */
-        void AddToPvs(const RayInfoContainer &rays, 
-                                  int &sampleContributions,
-                                  int &contributingSamples, 
-                                  bool storeLeavesWithRays = false);
-
-protected:
-
-        /// if NULL this does not correspond to feasible viewcell
-        BspViewCell *mViewCell;
-};
 
 /**
@@ -345 +44 @@
         {  
                 /// the current node
-                VspBspNode *mNode;
+                BspNode *mNode;
                 /// polygonal data for splitting
                 PolygonContainer *mPolygons;
                 /// current depth
                 int mDepth;
-                /// the view cell associated with this subdivsion
-                ViewCell *mViewCell;
+                
                 /// rays piercing this node
                 RayInfoContainer *mRays;
                 /// area of current node
                 float mArea;
-                /// geometry of current node induced by split planes
-                VspBspNodeGeometry *mGeometry;
-
+                /// geometry of node as induced by planes
+                BspNodeGeometry *mGeometry;
+                
                 /// pvs size
                 int mPvs;
@@ -374 +72 @@
                 mPolygons(NULL),
                 mDepth(0),
-                mViewCell(NULL),
                 mRays(NULL),
                 mPvs(0),
@@ -381 +78 @@
                 {}
                 
-                VspBspTraversalData(VspBspNode *node, 
+                VspBspTraversalData(BspNode *node, 
                                                  PolygonContainer *polys, 
                                                  const int depth, 
-                                                 ViewCell *viewCell,
                                                  RayInfoContainer *rays,
                                                  int pvs,
                                                  float area,
-                                                 VspBspNodeGeometry *cell): 
+                                                 BspNodeGeometry *geom): 
                 mNode(node), 
                 mPolygons(polys), 
                 mDepth(depth), 
-                mViewCell(viewCell),
                 mRays(rays),
                 mPvs(pvs),
                 mArea(area),
-                mGeometry(cell)
+                mGeometry(geom)
+                {}
+
+                VspBspTraversalData(PolygonContainer *polys, 
+                                                        const int depth, 
+                                                        RayInfoContainer *rays,
+                                                        BspNodeGeometry *geom): 
+                mNode(NULL), 
+                mPolygons(polys), 
+                mDepth(depth), 
+                mRays(rays),
+                mPvs(0),
+                mArea(0),
+                mGeometry(geom)
                 {}
     };
@@ -406 +114 @@
         VspBspTree();
 
+        /** Default destructor.
+        */
         ~VspBspTree();
 
-        const VspBspTreeStatistics &GetStatistics() const; 
+        /** Returns BSP Tree statistics.
+        */
+        const BspTreeStatistics &GetStatistics() const; 
   
 
@@ -420 +132 @@
         /** Returns list of BSP leaves.
         */
-        void CollectLeaves(vector<VspBspLeaf *> &leaves) const;
+        void CollectLeaves(vector<BspLeaf *> &leaves) const;
 
         /** Returns box which bounds the whole tree.
@@ -428 +140 @@
         /** Returns root of BSP tree.
         */
-        VspBspNode *GetRoot() const;
+        BspNode *GetRoot() const;
 
         /** Exports VspBsp tree to file.
@@ -450 +162 @@
         /** Returns statistics.
         */
-        VspBspTreeStatistics &GetStat();
+        BspTreeStatistics &GetStat();
 
         /** finds neighbouring leaves of this tree node.
         */
-        int FindNeighbors(VspBspNode *n, vector<VspBspLeaf *> &neighbors, 
+        int FindNeighbors(BspNode *n, 
+                                          vector<BspLeaf *> &neighbors, 
                                           const bool onlyUnmailed) const;
 
@@ -460 +173 @@
                 leading to this node.
         */
-        void ConstructGeometry(VspBspNode *n, PolygonContainer &cell) const;
+        void ConstructGeometry(BspNode *n, PolygonContainer &cell) const;
 
         /** Constructs geometry associated with the half space intersections 
@@ -469 +182 @@
         /** Construct geometry and stores it in a geometry node container.
         */
-        void ConstructGeometry(VspBspNode *n, VspBspNodeGeometry &cell) const;
+        void ConstructGeometry(BspNode *n, BspNodeGeometry &cell) const;
 
         /** Returns random leaf of BSP tree.
                 @param halfspace defines the halfspace from which the leaf is taken.
         */
-        VspBspLeaf *GetRandomLeaf(const Plane3 &halfspace);
+        BspLeaf *GetRandomLeaf(const Plane3 &halfspace);
 
         /** Returns random leaf of BSP tree.
                 @param onlyUnmailed if only unmailed leaves should be returned.
         */
-        VspBspLeaf *GetRandomLeaf(const bool onlyUnmailed = false);
+        BspLeaf *GetRandomLeaf(const bool onlyUnmailed = false);
 
         /** Traverses tree and counts all view cells as well as their PVS size.
         */
-        void EvaluateViewCellsStats(VspBspViewCellsStatistics &stat) const;
+        void EvaluateViewCellsStats(BspViewCellsStatistics &stat) const;
 
 
@@ -489 +202 @@
         */
         BspViewCell *GetRootCell() const;
+
+        /** Returns epsilon of this tree.
+        */
+        float GetEpsilon() const;
 
 protected:
@@ -521 +238 @@
                 @returns new root of the subtree
         */
-        VspBspNode *Subdivide(VspBspTraversalStack &tStack, VspBspTraversalData &tData);
-
-        /** Constructs the tree from the given list of polygons and rays.
+        BspNode *Subdivide(VspBspTraversalStack &tStack, 
+                                           VspBspTraversalData &tData);
+
+        /** Constructs the tree from the given traversal data.
                 @param polys stores set of polygons on which subdivision may be based
                 @param rays storesset of rays on which subdivision may be based
         */
-        void Construct(PolygonContainer *polys, RayInfoContainer *rays);
+        void Construct(const PolygonContainer &polys, RayInfoContainer *rays);
 
         /** Selects the best possible splitting plane. 
@@ -536 +254 @@
                 @returns the split plane
         */
-        Plane3 SelectPlane(VspBspLeaf *leaf, 
+        Plane3 SelectPlane(BspLeaf *leaf, 
                                            VspBspTraversalData &data);
 
@@ -564 +282 @@
         */
 
-        VspBspInterior *SubdivideNode(VspBspTraversalData &tData,
+        BspInterior *SubdivideNode(VspBspTraversalData &tData,
                                                            VspBspTraversalData &frontData,
                                                            VspBspTraversalData &backData,
@@ -575 +293 @@
                 @param rays bundle of rays on which the split can be based
         */
-        Plane3 SelectPlaneHeuristics(VspBspLeaf *leaf,
+        Plane3 SelectPlaneHeuristics(BspLeaf *leaf,
                                                                  VspBspTraversalData &data);
 
@@ -639 +357 @@
         bool SelectAxisAlignedPlane(Plane3 &plane, const PolygonContainer &polys) const;
 
-        /** Bounds ray and returns minT and maxT.
-                @returns true if ray hits BSP tree bounding box
-        */
-        bool BoundRay(const Ray &ray, float &minT, float &maxT) const;
-
         /** Subdivides the rays into front and back rays according to the split plane.
                 
@@ -662 +375 @@
         /** Extracts the split planes representing the space bounded by node n.
         */
-        void ExtractHalfSpaces(VspBspNode *n, vector<Plane3> &halfSpaces) const;
+        void ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const;
 
         /** Adds the object to the pvs of the front and back leaf with a given classification.
@@ -686 +399 @@
         float AccumulatedRayLength(const RayInfoContainer &rays) const;
 
-
+        /** Splits polygons with respect to the split plane.
+
+                @param plane the split plane
+                @param polys the polygons to be split. the polygons are consumed and
+                           distributed to the containers frontPolys, backPolys, coincident.
+                @param frontPolys returns the polygons in the front of the split plane
+                @param backPolys returns the polygons in the back of the split plane
+                @param coincident returns the polygons coincident to the split plane
+
+                @returns the number of splits   
+        */
+        int SplitPolygons(const Plane3 &plane,
+                                          PolygonContainer &polys, 
+                                          PolygonContainer &frontPolys, 
+                                          PolygonContainer &backPolys, 
+                                          PolygonContainer &coincident) const;
+
+        /** Adds ray sample contributions to the PVS.
+                @param sampleContributions the number contributions of the samples
+                @param contributingSampels the number of contributing rays
+                
+        */
+        void AddToPvs(BspLeaf *leaf,
+                                  const RayInfoContainer &rays, 
+                                  int &sampleContributions,
+                                  int &contributingSamples);
 
         /// Pointer to the root of the tree
-        VspBspNode *mRoot;
-                
-        VspBspTreeStatistics mStat;
+        BspNode *mRoot;
+                
+        BspTreeStatistics mStat;
 
         /// Strategies for choosing next split plane.
@@ -749 +487 @@
         bool mPvsUseArea;
 
+        float mEpsilon;
+
+        int mMaxTests;
+
 private:
         

trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp

@@ -459 +459 @@
         
         mStat.nodes = 1;
-        
-        mBox.Initialize();
-        
-        for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
-        {
-                leaf->AddRay(RayInfo(*ri));
-
-                mBox.Include((*ri)->GetOrigin());
-                mBox.Include((*ri)->GetTermination()); 
-        }
-        
+        mBox.Initialize();
+
+        //-- compute bounding box
         if (forcedBoundingBox)
                 mBox = *forcedBoundingBox;
+        else
+                for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
+                {
+                        mBox.Include((*ri)->GetOrigin());
+                        mBox.Include((*ri)->GetTermination()); 
+                }
 
         Debug << "bbox = " << mBox << endl;
+
+        for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
+        {
+                //leaf->AddRay(RayInfo(*ri));
+                VssRay *ray = *ri;
+                
+                float minT, maxT;
+
+                // TODO: not very efficient to implictly cast between rays types ...
+                if (mBox.GetRaySegment(*ray, minT, maxT))
+                {
+                        float len = ray->Length();
+                        
+                        if (!len) 
+                                len = Limits::Small;
+                
+                        leaf->AddRay(RayInfo(ray, minT / len, maxT / len)); 
+                }
+        }
         
         mStat.rays = (int)leaf->mRays.size();
@@ -1317 +1334 @@
                                                                          stack<RayTraversalData> &tstack)
 {
-        VspKdTreeInterior *in =  (VspKdTreeInterior *) data.mNode;
+        VspKdTreeInterior *in = dynamic_cast<VspKdTreeInterior *>(data.mNode);
   
         if (in->mAxis <= VspKdTreeNode::SPLIT_Z) 
@@ -1688 +1705 @@
     
         AxisAlignedBox3 box(node->mParent->mBox);
-        if (node->mParent->mFront == node)
+        if (node->mParent->GetFront() == node)
                 box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
         else
@@ -1763 +1780 @@
                         VspKdTreeInterior *interior = dynamic_cast<VspKdTreeInterior *>(node);
 
-                        nodeStack.push(interior->mBack);
-                        nodeStack.push(interior->mFront);
-                }
-        }
-}
+                        nodeStack.push(interior->GetBack());
+                        nodeStack.push(interior->GetFront());
+                }
+        }
+}

trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp

@@ -384 +384 @@
   VssTree *vssTree = NULL;
 
-  RayContainer bspRays;
-
   while (totalSamples < mInitialSamples) {
         int passContributingSamples = 0;
@@ -454 +452 @@
   mSceneGraph->CollectObjects(&objects);
 
-  //mViewCellsManager->Construct(mVssRays, objects);
+  // construct view cells
+  mViewCellsManager->Construct(objects, mVssRays, mViewSpaceBox);
 
   vssTree = new VssTree;
@@ -489 +488 @@
   int samples = 0;
   int pass = 0;
+
+  /// Rays used for post processing and visualizations.
+  RayContainer storedRays;
+  // cast view cell samples
   while (1) {
         int num = mVssSamplesPerPass;
@@ -503 +506 @@
           num = GenerateImportanceRays(vssTree, num, rays);
         }
-                
-                
+                        
         for (int i=0; i < rays.size(); i++)
           CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays);
@@ -528 +530 @@
         }
 
-        // cast rays into BSP tree
-        /*if (ViewCell::sHierarchy == ViewCell::BSP) 
-          {
-                for (int i = 0; i < (int)vssRays.size(); ++ i)
-                  {
-                        CastRay(*mBspTree, *vssRays[i]);
-                  }
-          }
-*/
+        //-- prepare traversal rays for view cell intersections
+        RayContainer passRays;
+
+        VssRayContainer::const_iterator it, it_end = vssRays.end();
+        
+        for (it = vssRays.begin(); it != it_end; ++ it)
+                passRays.push_back(new Ray(*(*it)));
+        
+        int sampleContributions = 0;
+        int contributingSamples = 0;
+
+        /// compute view cell contribution of rays
+        mViewCellsManager->ComputeSampleContributions(passRays,
+                                                                                                  sampleContributions, 
+                                                                                                  contributingSamples);
+        
+        //-- save rays for post processing
+        if (((int)storedRays.size() < mViewCellsManager->GetPostProcessSamples()) ||
+                ((int)storedRays.size() < mViewCellsManager->GetVisualizationSamples()))
+        {
+                RayContainer::const_iterator it, it_end = passRays.end();
+
+                for (it = passRays.begin(); it != it_end; ++ it)
+                        storedRays.push_back(new Ray(*(*it)));
+        }
+        else
+        {
+                CLEAR_CONTAINER(passRays);
+        }
+
         samples+=num;
         float pvs = vssTree->GetAvgPvsSize();
@@ -547 +570 @@
         pass++;
   }
-        
+
+  //-- post process view cells
+  mViewCellsManager->PostProcess(objects, storedRays);
+  mViewCellsManager->Visualize(objects, storedRays);
+
+  CLEAR_CONTAINER(storedRays);
+
   delete vssTree;
 

trunk/VUT/GtpVisibilityPreprocessor/src/VssRay.h

@@ -135 +135 @@
   Vector3 GetNormalizedDir() const { return (mTermination - mOrigin)*mInvSize; }
 
+  float Length() const { return Distance(mOrigin, mTermination); }
+
   Vector3 Extrap(const float t) const {
         return GetOrigin() + t * GetDir();

trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h

@@ -19 +19 @@
 class ViewCell;
 class BspTree;
+class VspBspTree;
+class BspNode;
 
 class X3dExporter : public Exporter
@@ -92 +94 @@
   ExportGeometry(const ObjectContainer &objects);
 
+  bool
+  ExportRays(const RayContainer &rays,
+                         const float length=1000,
+                         const RgbColor &color = RgbColor(1,1,1));
+  bool
+  ExportRays(const VssRayContainer &rays,
+                         const RgbColor &color = RgbColor(1,1,1));
+
   virtual void 
   ExportBspSplitPlanes(const BspTree &tree);
@@ -101 +111 @@
   ExportLeavesGeometry(const BspTree &tree, const vector<BspLeaf *> &leaves);
 
-  bool
-  ExportRays(const RayContainer &rays,
-                                                 const float length=1000,
-                                                 const RgbColor &color = RgbColor(1,1,1));
-
-        bool
-        ExportRays(const VssRayContainer &rays,
-                                                 const RgbColor &color = RgbColor(1,1,1));
-                
   virtual void
   ExportBspViewCellPartition(const BspTree &tree, const int maxPvs = 0);
@@ -116 +117 @@
   ExportBspLeaves(const BspTree &tree, const int maxPvs = 0);
 
+
+  virtual void 
+  ExportBspSplits(const VspBspTree &tree, const bool exportDepth);
+
+  virtual void
+  ExportBspViewCellPartition(const VspBspTree &tree, const int maxPvs = 0);
+
 protected:
+
 
   virtual void
@@ -126 +135 @@
   bool
   ExportBspTreeRayDensity(const BspTree &tree);
+  
+  void
+  AddBoxToMesh(const AxisAlignedBox3 &box,
+                           Mesh *mesh);
 
-        void
-        AddBoxToMesh(const AxisAlignedBox3 &box,
-                                                         Mesh *mesh);
-                
+  void ExportBspNodeSplits(BspNode *root,
+                                                   const AxisAlignedBox3 &box, 
+                                                   const bool exportDepth,
+                                                   const bool epsilon);
 };
  

trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp

@@ -58 +58 @@
         p->KdTreeStatistics(cout);
         
-        // parse view cell hierarchy options
-        p->ParseViewCellsOptions();
-  
         // parse view cells related options
         p->PrepareViewCells();