Changeset 508 for trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp

Timestamp:

01/08/06 05:56:40 (19 years ago)

Author:

mattausch

Message:

implemented view cells exporting / loading
improved vsp bsp tree (only axis aligbed until a level), reuse results from Plane
testing, collectmergeneighbors
implemented view cell meshes

File:

• trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp (modified) (45 diffs)

trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp

@@ -19 +19 @@
 
 //-- static members
+
 /** Evaluates split plane classification with respect to the plane's
         contribution for a minimum number of ray splits.
@@ -35 +36 @@
 float BspMergeCandidate::sOverallCost = 0;
 
-/****************************************************************/
-/*                  class VspBspTree implementation             */
-/****************************************************************/
+/********************************************************************/
+/*                  class VspBspTree implementation                 */
+/********************************************************************/
 
 VspBspTree::VspBspTree():
@@ -44 +45 @@
 mCostNormalizer(Limits::Small),
 mViewCellsManager(NULL),
-mStoreRays(false),
 mOutOfBoundsCell(NULL),
-mShowInvalidSpace(false)
+mShowInvalidSpace(false),
+mStoreRays(false)
 {
         bool randomize = false;
@@ -88 +89 @@
         environment->GetIntValue("ViewCells.maxPvs", mMaxPvs);
 
+        //-- termination criteria for axis aligned split
+        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxRayContribution", 
+                mTermMaxRayContriForAxisAligned);
+        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
+                mTermMinRaysForAxisAligned);
+
+        //environment->GetFloatValue("VspBspTree.maxTotalMemory", mMaxTotalMemory);
+        environment->GetFloatValue("VspBspTree.maxStaticMemory", mMaxMemory);
 
         //-- debug output
@@ -135 +144 @@
 }
 
+BspViewCell *VspBspTree::GetOutOfBoundsCell()
+{
+        return mOutOfBoundsCell;
+}
+
 
 BspViewCell *VspBspTree::GetOrCreateOutOfBoundsCell()
 {
         if (!mOutOfBoundsCell)
+        {
                 mOutOfBoundsCell = new BspViewCell();
-
+                mOutOfBoundsCell->SetId(-1);
+        }
         return mOutOfBoundsCell;
 }
@@ -326 +342 @@
 }
 
+
+// return memory usage in MB
+float VspBspTree::GetMemUsage() const
+{
+        return
+                (sizeof(VspBspTree) +
+                 mStat.Leaves() * sizeof(BspLeaf) +
+                 mStat.Interior() * sizeof(BspInterior) +
+                 mStat.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
+}
+
+
+
 void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
 {
@@ -350 +379 @@
 
         long startTime = GetTime();
+        
+        bool mOutOfMemory = false;
 
         while (!tStack.empty())
         {
                 tData = tStack.top();
-
-            tStack.pop();
+            tStack.pop();               
+
+                if (0 && !mOutOfMemory)
+                {
+                        float mem = GetMemUsage();
+
+                        if (mem > mMaxMemory)
+                        {
+                                mOutOfMemory = true;
+                                Debug << "memory limit reached: " << mem << endl;
+                        }
+                }
 
                 // subdivide leaf node
@@ -365 +406 @@
         }
 
+        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl;
         cout << "finished\n";
 
         mStat.Stop();
 }
+
 
 bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
@@ -377 +420 @@
                  (data.mArea <= mTermMinArea) ||
                  (mStat.Leaves() >= mMaxViewCells) ||
-                // (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
+                 (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
                  (data.mDepth >= mTermMaxDepth));
 }
+
 
 BspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack,
@@ -386 +430 @@
         BspNode *newNode = tData.mNode;
 
-        if (!TerminationCriteriaMet(tData))
+        if (!mOutOfMemory || !TerminationCriteriaMet(tData))
         {
                 PolygonContainer coincident;
@@ -423 +467 @@
                 }
                 else
-                {
-                        // create new view cell for this leaf
+                {       // create new view cell for this leaf                   
                         viewCell = new BspViewCell();
                 }
@@ -465 +508 @@
 {
         BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
-
-        int maxCostMisses = tData.mMaxCostMisses;
-
+        
         // select subdivision plane
         Plane3 splitPlane;
-        if (!SelectPlane(splitPlane, leaf, tData))
+        int maxCostMisses = tData.mMaxCostMisses;
+        bool success = SelectPlane(splitPlane, 
+                                                           leaf, 
+                                                           tData,
+                                                           frontData,
+                                                           backData);
+        if (!success)
         {
                 ++ maxCostMisses;
@@ -492 +539 @@
 
         //-- the front and back traversal data is filled with the new values
+        frontData.mDepth = tData.mDepth + 1;
         frontData.mPolygons = new PolygonContainer();
-        frontData.mDepth = tData.mDepth + 1;
         frontData.mRays = new RayInfoContainer();
-        frontData.mGeometry = new BspNodeGeometry();
-
+        
+        backData.mDepth = tData.mDepth + 1;
         backData.mPolygons = new PolygonContainer();
-        backData.mDepth = tData.mDepth + 1;
         backData.mRays = new RayInfoContainer();
-        backData.mGeometry = new BspNodeGeometry();
-
+        
         // subdivide rays
         SplitRays(interior->GetPlane(),
@@ -524 +569 @@
         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);
+        // split front and back node geometry and compute area
+        if (mPvsUseArea)
+        {
+                // if not already computed
+                if (!frontData.mGeometry && !backData.mGeometry)
+                {
+                        frontData.mGeometry = new BspNodeGeometry();
+                        backData.mGeometry = new BspNodeGeometry();
+
+                        tData.mGeometry->SplitGeometry(*frontData.mGeometry,
+                                                                                   *backData.mGeometry,
+                                                                                   interior->GetPlane(),
+                                                                                   mBox,
+                                                                                   mEpsilon);
+                
+                        frontData.mArea = frontData.mGeometry->GetArea();
+                        backData.mArea = backData.mGeometry->GetArea();
+                }
+        }
+        else
+        {
+                frontData.mArea = (float)frontData.mRays->size();
+                backData.mArea = (float)backData.mRays->size();
+        }
 
         //-- create front and back leaf
@@ -565 +618 @@
         return interior;
 }
+
 
 void VspBspTree::AddToPvs(BspLeaf *leaf,
@@ -757 +811 @@
                                                                                  const VspBspTraversalData &tData,
                                                                                  int &axis,
-                                                                                 float &position,
-                                                                                 int &raysBack,
-                                                                                 int &raysFront,
-                                                                                 int &pvsBack,
-                                                                                 int &pvsFront)
-{
-        AxisAlignedBox3 box;
-        box.Initialize();
-
-        // create bounding box of region
-        RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();
-
-        for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
-                box.Include((*ri).ExtrapTermination());
-
+                                                                                 BspNodeGeometry **frontGeom,
+                                                                                 BspNodeGeometry **backGeom,
+                                                                                 float &frontArea,
+                                                                                 float &backArea)
+{
         const bool useCostHeuristics = false;
 
@@ -777 +821 @@
         float nPosition[3];
         float nCostRatio[3];
+        float nFrontArea[3];
+        float nBackArea[3];
+
+        BspNodeGeometry *nFrontGeom[3];
+        BspNodeGeometry *nBackGeom[3];
+
         int bestAxis = -1;
 
-        const int sAxis = box.Size().DrivingAxis();
+        // create bounding box of node extent
+        AxisAlignedBox3 box;
+        box.Initialize();
+        
+#if 0
+        RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();
+
+        for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
+                box.Include((*ri).ExtrapTermination());
+#else
+        PolygonContainer::const_iterator it, it_end = tData.mGeometry->mPolys.end();
+
+        for(it = tData.mGeometry->mPolys.begin(); it < it_end; ++ it)
+                box.Include(*(*it));
+
+#endif
+        int sAxis = box.Size().DrivingAxis();
 
         for (axis = 0; axis < 3; ++ axis)
         {
+                nFrontGeom[axis] = new BspNodeGeometry();
+                nBackGeom[axis] = new BspNodeGeometry();
+
                 if (!mOnlyDrivingAxis || axis == sAxis)
                 {
@@ -789 +858 @@
                                 nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
                                 Vector3 normal(0,0,0); normal[axis] = 1;
-
-                                nCostRatio[axis] = SplitPlaneCost(Plane3(normal, nPosition[axis]), tData);
+                                
+                                nCostRatio[axis] = SplitPlaneCost(Plane3(normal, nPosition[axis]), 
+                                                                                                  tData, *nFrontGeom[axis], *nBackGeom[axis],
+                                                                                                  nFrontArea[axis], nBackArea[axis]);
                         }
                         else
@@ -809 +880 @@
                         else if (nCostRatio[axis] < nCostRatio[bestAxis])
                         {
-                                /*Debug << "pvs front " << nPvsBack[axis]
-                                          << " pvs back " << nPvsFront[axis]
-                                          << " overall pvs " << leaf->GetPvsSize() << endl;*/
                                 bestAxis = axis;
                         }
@@ -820 +888 @@
         //-- assign values
         axis = bestAxis;
-        position = nPosition[bestAxis];
-
-        raysBack = 0;//nRaysBack[bestAxis];
-        raysFront = 0;//nRaysFront[bestAxis];
-
-        pvsBack = 0;//nPvsBack[bestAxis];
-        pvsFront = 0;//nPvsFront[bestAxis];
-
-        Vector3 normal(0,0,0); normal[bestAxis] = 1;
+        frontArea = nFrontArea[bestAxis];
+        backArea = nBackArea[bestAxis];
+
+        // assign best split nodes geometry 
+        *frontGeom = nFrontGeom[bestAxis];
+        *backGeom = nBackGeom[bestAxis];
+        // and delete other geometry
+        delete nFrontGeom[(bestAxis + 1) % 3];
+        delete nBackGeom[(bestAxis + 2) % 3];
+
+        //-- split plane
+    Vector3 normal(0,0,0); normal[bestAxis] = 1;
         plane = Plane3(normal, nPosition[bestAxis]);
-        
+
         return nCostRatio[bestAxis];
 }
 
 
-/*
-float VspBspTree::EvalCostRatio(const VspBspTraversalData &tData,
-                                                                const AxisAlignedBox3 &box,
-                                                                const int axis,
-                                                                const float position,
-                                                                int &raysBack,
-                                                                int &raysFront,
-                                                                int &pvsBack,
-                                                                int &pvsFront)
-{
-        raysBack = 0;
-        raysFront = 0;
-        pvsFront = 0;
-        pvsBack = 0;
-
-        Intersectable::NewMail(3);
-
-        // eval pvs size
-        const int pvsSize = tData.mPvs;
-
-        // create unique ids for pvs heuristics
-        GenerateUniqueIdsForPvs();
-
-        // this is the main ray classification loop!
-        for(RayInfoContainer::iterator ri = tData.mRays->begin();
-                        ri != tData.mRays->end(); ++ ri)
-        {
-                if (!(*ri).mRay->IsActive())
-                        continue;
-
-                // determine the side of this ray with respect to the plane
-                float t;
-                int side = (*ri).ComputeRayIntersection(axis, position, t);
-                        
-                if (side <= 0)
-                        ++ raysBack;
-
-                if (side >= 0)
-                        ++ raysFront;
-
-                AddObjToPvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
-        }
-
-        //-- only one of these options should be one
-
-        if (0) //-- only pvs
-        {
-                const float sum = float(pvsBack + pvsFront);
-                const float oldCost = (float)pvsSize;
-
-                return sum / oldCost;
-        }
-
-        //-- pvs + probability heuristics
-        float pBack, pFront, pOverall;
-
-        if (0)
-        {
-                // box length substitute for probability
-                const float minBox = box.Min(axis);
-                const float maxBox = box.Max(axis);
-
-                pBack = position - minBox;
-                pFront = maxBox - position;
-                pOverall = maxBox - minBox;
-        }
-
-        if (1) //-- area substitute for probability
-        {
-                
-                const bool useMidSplit = true;
-                //const bool useMidSplit = false;
-                        
-                pOverall = box.SurfaceArea();
-                        
-                if (!useMidSplit)
-                {
-                        Vector3 pos = box.Max(); pos[axis] = position;
-                        pBack = AxisAlignedBox3(box.Min(), pos).SurfaceArea();
-
-                        pos = box.Min(); pos[axis] = position;
-                        pFront = AxisAlignedBox3(pos, box.Max()).SurfaceArea();
-                }
-                else
-                {
-                        //-- simplified computation for mid split
-                        const int axis2 = (axis + 1) % 3;
-                        const int axis3 = (axis + 2) % 3;
-
-                        const float faceArea = 
-                                (box.Max(axis2) - box.Min(axis2)) *
-                                (box.Max(axis3) - box.Min(axis3));
-
-                        pBack = pFront = pOverall * 0.5f + faceArea;
-                }
-        }
-
-        //-- ray density substitute for probability
-        if (0)
-        {
-                pBack = (float)raysBack;
-                pFront = (float)raysFront;
-                pOverall = (float)tData.mRays->size();
-        }
-
-        //Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
-        //Debug << pFront << " " << pBack << " " << pOverall << endl;
-
-        // float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
-        const float newCost = pvsBack * pBack + pvsFront * pFront;
-        //  float oldCost = leaf->mRays.size();
-        const float oldCost = (float)pvsSize * pOverall;
-
-        return  (mCtDivCi + newCost) / oldCost;
-}
-*/
-
-
-bool VspBspTree::SelectPlane(Plane3 &plane,
+bool VspBspTree::SelectPlane(Plane3 &bestPlane,
                                                          BspLeaf *leaf,
-                                                         VspBspTraversalData &data)
+                                                         VspBspTraversalData &data,
+                                                         VspBspTraversalData &frontData,
+                                                         VspBspTraversalData &backData)
 {
         // simplest strategy: just take next polygon
@@ -962 +917 @@
         if (!data.mPolygons->empty())
                 {
-                        const int randIdx = (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
+                        const int randIdx = 
+                                (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
                         Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
 
-                        plane = nextPoly->GetSupportingPlane();
+                        bestPlane = nextPoly->GetSupportingPlane();
                         return true;
                 }
         }
 
-        // use heuristics to find appropriate plane
-        return SelectPlaneHeuristics(plane, leaf, data);
+        //-- use heuristics to find appropriate plane
+
+        // intermediate plane
+        Plane3 plane;
+        float lowestCost = MAX_FLOAT;
+                
+        const bool onlyAxisAligned  = 
+                (mSplitPlaneStrategy & AXIS_ALIGNED) &&
+                ((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
+                ((int)data.GetAvgRayContribution() < mTermMaxRayContriForAxisAligned);
+        
+        // split polygons if no axis aligned splits
+        const int limit = onlyAxisAligned ? 0 : 
+                Min((int)data.mPolygons->size(), mMaxPolyCandidates);
+
+        float candidateCost;
+
+        int maxIdx = (int)data.mPolygons->size();
+
+        for (int i = 0; i < limit; ++ i)
+        {
+                //-- assure that no index is taken twice
+                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
+                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
+
+                // swap candidate to the end to avoid testing same plane
+                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
+                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
+
+                // evaluate current candidate
+                BspNodeGeometry fGeom, bGeom;
+                float fArea, bArea;
+                plane = poly->GetSupportingPlane();
+                candidateCost = SplitPlaneCost(plane, data, fGeom, bGeom, fArea, bArea);
+                
+                if (candidateCost < lowestCost)
+                {
+                        bestPlane = plane;
+                        lowestCost = candidateCost;
+                }
+        }
+
+        //-- evaluate axis aligned splits
+        int axis;
+        BspNodeGeometry *fGeom, *bGeom;
+        float fArea, bArea;
+
+        candidateCost = SelectAxisAlignedPlane(plane, data, axis,
+                                                                                   &fGeom, &bGeom, 
+                                                                                   fArea, bArea);                                                
+
+        if (candidateCost < lowestCost)
+        {       
+                bestPlane = plane;
+                lowestCost = candidateCost;
+
+                //-- assign already computed values
+                frontData.mGeometry = fGeom;
+                backData.mGeometry = bGeom;
+                frontData.mArea = fArea;
+                backData.mArea = bArea;
+
+                //! error also computed if cost ratio is missed
+                ++ mStat.splits[axis];
+        }
+        else
+        {
+                DEL_PTR(fGeom);
+                DEL_PTR(bGeom);
+        }
+
+#ifdef _DEBUG
+        Debug << "plane lowest cost: " << lowestCost << endl;
+#endif
+
+        // cost ratio miss
+        if (lowestCost > mTermMaxCostRatio)
+                return false;
+
+        return true;
 }
 
@@ -1048 +1082 @@
 
 
-bool VspBspTree::SelectPlaneHeuristics(Plane3 &bestPlane,
-                                                                           BspLeaf *leaf,
-                                                                           VspBspTraversalData &data)
-{
-        float lowestCost = MAX_FLOAT;
-        // intermediate plane
-        Plane3 plane;
-        bool useAxisAlignedPlane = false;
-
-        const int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
-        int maxIdx = (int)data.mPolygons->size();
-
-        float candidateCost;
-
-        for (int i = 0; i < limit; ++ i)
-        {
-                //-- assure that no index is taken twice
-                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
-                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
-
-                // swap candidate to the end to avoid testing same plane
-                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
-
-                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
-
-                // evaluate current candidate
-                candidateCost = SplitPlaneCost(poly->GetSupportingPlane(), data);
-
-                if (candidateCost < lowestCost)
-                {
-                        bestPlane = poly->GetSupportingPlane();
-                        lowestCost = candidateCost;
-                }
-        }
-
-        //-- axis aligned splits
-        int axis;
-        float position;
-        int raysBack;
-        int raysFront;
-        int pvsFront;
-        int pvsBack;
-
-        candidateCost = SelectAxisAlignedPlane(plane, 
-                                                                                   data, 
-                                                                                   axis,
-                                                                                   position,
-                                                                                   raysBack,
-                                                                                   raysFront,
-                                                                                   pvsFront,
-                                                                                   pvsBack);
-
-        if (candidateCost < lowestCost)
-        {       
-                if (!useAxisAlignedPlane)
-                {
-                        useAxisAlignedPlane = true;
-                        //if (data.mPolygons->size() > 0)
-                        //      Debug << "haha" << endl;
-                        //! error also computed if cost ratio is missed
-                        ++ mStat.splits[axis];
-                }
-
-                bestPlane = plane;
-                lowestCost = candidateCost;
-        }
-
-#ifdef _DEBUG
-        Debug << "plane lowest cost: " << lowestCost << endl;
-#endif
-
-        // cost ratio miss
-        if (lowestCost > mTermMaxCostRatio)
-                return false;
-
-        return true;
-}
-
-
 inline void VspBspTree::GenerateUniqueIdsForPvs()
 {
@@ -1136 +1091 @@
 
 float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane,
-                                                                 const VspBspTraversalData &data) const
+                                                                 const VspBspTraversalData &data,
+                                                                 BspNodeGeometry &geomFront,
+                                                                 BspNodeGeometry &geomBack,
+                                                                 float &areaFront,
+                                                                 float &areaBack) const
 {
         float cost = 0;
@@ -1143 +1102 @@
         float sumRaySplits = 0;
 
-        int frontPvs = 0;
-        int backPvs = 0;
-
-        // probability that view point lies in child
+        int pvsFront = 0;
+        int pvsBack = 0;
+
+        // probability that view point lies in back / front node
         float pOverall = 0;
         float pFront = 0;
         float pBack = 0;
 
+        int raysFront = 0;
+        int raysBack = 0;
+        int totalPvs = 0;
+
         if (mSplitPlaneStrategy & PVS)
         {
@@ -1159 +1122 @@
                 {
                         // construct child geometry with regard to the candidate split plane
-                        BspNodeGeometry frontCell;
-                        BspNodeGeometry backCell;
-
-                        data.mGeometry->SplitGeometry(frontCell,
-                                                                                  backCell,
+                        data.mGeometry->SplitGeometry(geomFront,
+                                                                                  geomBack,
                                                                                   candidatePlane,
                                                                                   mBox,
                                                                                   mEpsilon);
 
-                        pFront = frontCell.GetArea();
-                        pBack = backCell.GetArea();
-
+                        areaFront = geomFront.GetArea();
+                        areaBack = geomBack.GetArea();
+
+                        pBack = areaBack;
+                        pFront = areaFront;
                         pOverall = data.mArea;
+                }
+                else // use number of rays to approximate volume
+                {
+                        pOverall = (float)data.mRays->size();
+                        pFront = (float)raysFront;
+                        pBack = (float)raysBack;
                 }
         }
@@ -1189 +1157 @@
                 limit = (int)data.mRays->size();
         }
-
+        
         for (int i = 0; i < limit; ++ i)
         {
-                const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;
+                const int testIdx = useRand ? 
+                        (int)RandomValue(0, (Real)((int)data.mRays->size() - 1)) : i;
                 RayInfo rayInf = (*data.mRays)[testIdx];
 
+                float t;
                 VssRay *ray = rayInf.mRay;
-                float t;
                 const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
+
+        if (cf >= 0)
+                        ++ raysFront;
+                if (cf <= 0)
+                        ++ raysBack;
 
                 if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
@@ -1212 +1186 @@
                 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
-                        AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
-                        AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
-
-                        // use number of rays to approximate volume
-                        if (!mPvsUseArea)
-                        {
-                                pOverall = (float)data.mRays->size();
-
-                                if (cf >= 0)
-                                        ++ pFront;
-                                if (cf <= 0)
-                                        ++ pBack;
-                        }
+                        // find front and back pvs for origing and termination object
+                        AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
+                        AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
                 }
         }
@@ -1242 +1203 @@
         if (mSplitPlaneStrategy & PVS)
         {
-                const float oldCost = pOverall * (float)data.mPvs + Limits::Small;
-                cost += mPvsFactor * (frontPvs * pFront + backPvs * pBack) / oldCost;
-
-                //cost += mPvsFactor * 0.5 * (frontPvs * pFront + backPvs * pBack) / oldCost;
-                //Debug << "new cost: " << cost << " over" << frontPvs * pFront + backPvs * pBack << " old cost " << oldCost << endl;
-
-                if (0) // give penalty to unbalanced split
-                        if (((pFront * 0.2 + Limits::Small) > pBack) ||
-                                (pFront < (pBack * 0.2 + Limits::Small)))
-                                        cost += 0.5;
+                const float oldCost = pOverall * (float)totalPvs + Limits::Small;
+                cost += mPvsFactor * (pvsFront * pFront + pvsBack * pBack) / oldCost;
         }
 
 #ifdef _DEBUG
         Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
-                  << " frontpvs: " << frontPvs << " pFront: " << pFront
-                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
+                  << " frontpvs: " << pvsFront << " pFront: " << pFront
+                  << " backpvs: " << pvsBack << " pBack: " << pBack << endl << endl;
 #endif
 
@@ -1263 +1216 @@
         return cost / (float)mCostNormalizer;
 }
+
 
 void VspBspTree::AddObjToPvs(Intersectable *obj,
                                                          const int cf,
                                                          int &frontPvs,
-                                                         int &backPvs) const
+                                                         int &backPvs,
+                                                         int &totalPvs) const
 {
         if (!obj)
                 return;
+        
+        if ((obj->mMailbox != sFrontId) &&
+                (obj->mMailbox != sBackId) &&
+                (obj->mMailbox != sFrontAndBackId))
+        {
+                ++ totalPvs;
+        }
+
         // TODO: does this really belong to no pvs?
         //if (cf == Ray::COINCIDENT) return;
@@ -1281 +1244 @@
                 {
                         ++ frontPvs;
-
+                
                         if (obj->mMailbox == sBackId)
                                 obj->mMailbox = sFrontAndBackId;
@@ -1295 +1258 @@
                 {
                         ++ backPvs;
-
+                
                         if (obj->mMailbox == sFrontId)
                                 obj->mMailbox = sFrontAndBackId;
@@ -1322 +1285 @@
                         BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
                         if (leaf->TreeValid() && 
-                                (!onlyUnmailed || leaf->Mailed()) &&
+                                (!onlyUnmailed || !leaf->Mailed()) &&
                                 ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
                         {
@@ -1351 +1314 @@
 
 
-BspViewCell *VspBspTree::GetOutOfBoundsCell() const
-{
-        return mOutOfBoundsCell;
-}
-
-
 void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
 {
@@ -1368 +1325 @@
         if (data.mPvs > mStat.maxPvs)
                 mStat.maxPvs = data.mPvs;
+        
         if (data.mDepth < mStat.minDepth)
                 mStat.minDepth = data.mDepth;
+
         if (data.mDepth >= mTermMaxDepth)
                 ++ mStat.maxDepthNodes;
+        // accumulate rays to compute rays /  leaf
+        mStat.accumRays += (int)data.mRays->size();
 
         if (data.mPvs < mTermMinPvs)
@@ -1383 +1344 @@
 
         if (data.mArea <= mTermMinArea)
-        {
-                //Debug << "area: " << data.mArea / mBox.SurfaceArea() << " min area: " << mTermMinArea / mBox.SurfaceArea() << endl;
                 ++ mStat.minAreaNodes;
-        }
+        
         // accumulate depth to compute average depth
         mStat.accumDepth += data.mDepth;
@@ -1498 +1457 @@
 }
 
-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
 {
@@ -1516 +1462 @@
 
         if (!mRoot)
-        return;
+                return;
 
         nodeStack.push(mRoot);
@@ -1543 +1489 @@
                 {
                         BspInterior *interior = dynamic_cast<BspInterior *>(node);
-                        
+                
                         nodeStack.push(interior->GetFront());
                         nodeStack.push(interior->GetBack());
@@ -2245 +2191 @@
 {
         BspLeaf::NewMail();
-
+        
         vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
 
+        Debug << "mergequeue size: " << mMergeQueue.size() << endl;
         // find merge candidates and push them into queue
         for (it = leaves.begin(); it != it_end; ++ it)
@@ -2268 +2215 @@
                 // TODO: test if at least one ray goes from one leaf to the other
                 for (nit = neighbors.begin(); nit != nit_end; ++ nit)
-                {                       
-                        mMergeQueue.push(BspMergeCandidate(leaf, *nit));
-                }
-        }
+                {
+                        if ((*nit)->GetViewCell() != leaf->GetViewCell())
+                                mMergeQueue.push(BspMergeCandidate(leaf, *nit));
+                }
+        }
+
+        Debug << "new mergequeue size: " << mMergeQueue.size() << endl;
 
         return (int)leaves.size();
@@ -2283 +2233 @@
         long startTime = GetTime();
         ConstructBspRays(bspRays, rays);
-        Debug << (int)bspRays.size() << " bsp rays constructed in " << TimeDiff(startTime, GetTime()) * 1e-3f << " secs" << endl;
+        Debug << (int)bspRays.size() << " bsp rays constructed in " 
+                  << TimeDiff(startTime, GetTime()) * 1e-3f << " secs" << endl;
 
         map<BspLeaf *, vector<BspLeaf*> > neighborMap;
@@ -2324 +2275 @@
             leaf = (*iit).mLeaf;
 
-                        // view space not valid
-                        if (!leaf->TreeValid() || !prevLeaf->TreeValid())
+                        // view space not valid or same view cell
+                        if (!leaf->TreeValid() || !prevLeaf->TreeValid() ||
+                                (leaf->GetViewCell() == prevLeaf->GetViewCell()))
                                 continue;
 
@@ -2379 +2331 @@
         vector<BspLeaf *> leaves;
         CollectLeaves(leaves, true, mMaxPvs);
-        Debug << "found " << (int)leaves.size() << " new leaves" << endl;
-        //CollectMergeCandidates(leaves);
+        Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
+        // TODO
+        CollectMergeCandidates(leaves);
 
         return numLeaves;
@@ -2508 +2461 @@
 
 
-ViewCell *
-VspBspTree::GetViewCell(const Vector3 &point)
+ViewCell *VspBspTree::GetViewCell(const Vector3 &point)
 {
   if (mRoot == NULL)
@@ -2737 +2689 @@
 }
 
-#define USE_ASCII 0
-bool VspBspTree::WriteVspBspTree()
-{
-        char fileName[100];
-        environment->GetStringValue("VspBspTree.viewCellsFilename", fileName);
+
+bool VspBspTree::Export(ofstream &stream)
+{
+        ExportNode(mRoot, stream);
+
+        return true;
+}
+
+
+void VspBspTree::ExportNode(BspNode *node, ofstream &stream)
+{
+        if (node->IsLeaf())
+        {
+                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
+                        
+                int id = -1;
+                if (leaf->GetViewCell() != mOutOfBoundsCell)
+                        id = leaf->GetViewCell()->GetId();
+
+                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
+        }
+        else
+        {
+                BspInterior *interior = dynamic_cast<BspInterior *>(node);
         
-        /*VssRayContainer::const_iterator it, it_end = samples.end();
-        
-#if USE_ASCII
-        ofstream samplesOut(fileName);
-        if (!samplesOut.is_open())
-                return false;
-
-        for (it = samples.begin(); it != it_end; ++ it)
-        {
-                VssRay *ray = *it;
-                int sourceid = ray->mOriginObject ? ray->mOriginObject->mId : -1;               
-                int termid = ray->mTerminationObject ? ray->mTerminationObject->mId : -1;       
-
-                samplesOut << ray->GetOrigin().x << " " << ray->GetOrigin().y << " " << ray->GetOrigin().z << " " 
-                                   << ray->GetTermination().x << " " << ray->GetTermination().y << " " << ray->GetTermination().z << " " 
-                                   << sourceid << " " << termid << "\n";
-        }
-#else
-        ofstream samplesOut(fileName, ios::binary);
-        if (!samplesOut.is_open())
-                return false;
-
-        for (it = samples.begin(); it != it_end; ++ it)
-        {       
-                VssRay *ray = *it;
-                Vector3 origin(ray->GetOrigin());
-                Vector3 termination(ray->GetTermination());
-                
-                int sourceid = ray->mOriginObject ? ray->mOriginObject->mId : -1;               
-                int termid = ray->mTerminationObject ? ray->mTerminationObject->mId : -1;               
-
-                samplesOut.write(reinterpret_cast<char *>(&origin), sizeof(Vector3));
-                samplesOut.write(reinterpret_cast<char *>(&termination), sizeof(Vector3));
-                samplesOut.write(reinterpret_cast<char *>(&sourceid), sizeof(int));
-                samplesOut.write(reinterpret_cast<char *>(&termid), sizeof(int));
-    }
-#endif
-
-        samplesOut.close();
-        */
-        return true;
-}
-
-bool VspBspTree::LoadVspBspTree()
-{
-        /*std::stable_sort(objects.begin(), objects.end(), ilt);
-        char fileName[100];
-        environment->GetStringValue("Preprocessor.samplesFilename", fileName);
-        
-    Vector3 origin, termination;
-        // HACK: needed only for lower_bound algorithm to find the 
-        // intersected objects
-        MeshInstance sObj(NULL);
-        MeshInstance tObj(NULL);
-
-#if USE_ASCII
-        ifstream samplesIn(fileName, ios::binary);
-        if (!samplesIn.is_open())
-                return false;
-
-        string buf;
-        while (!(getline(samplesIn, buf)).eof())
-        {
-                sscanf(buf.c_str(), "%f %f %f %f %f %f %d %d", 
-                           &origin.x, &origin.y, &origin.z,
-                           &termination.x, &termination.y, &termination.z, 
-                           &(sObj.mId), &(tObj.mId));
-                
-                Intersectable *sourceObj = NULL;
-                Intersectable *termObj = NULL;
-                
-                if (sObj.mId >= 0)
-                {
-                        ObjectContainer::iterator oit =
-                                lower_bound(objects.begin(), objects.end(), &sObj, ilt);
-                        sourceObj = *oit;
-                }
-                
-                if (tObj.mId >= 0)
-                {
-                        ObjectContainer::iterator oit =
-                                lower_bound(objects.begin(), objects.end(), &tObj, ilt);
-                        termObj = *oit;
-                }
-
-                samples.push_back(new VssRay(origin, termination, sourceObj, termObj));
-        }
-#else
-        ifstream samplesIn(fileName, ios::binary);
-        if (!samplesIn.is_open())
-                return false;
-
-        while (1)
-        {
-                 samplesIn.read(reinterpret_cast<char *>(&origin), sizeof(Vector3));
-                 samplesIn.read(reinterpret_cast<char *>(&termination), sizeof(Vector3));
-                 samplesIn.read(reinterpret_cast<char *>(&(sObj.mId)), sizeof(int));
-                 samplesIn.read(reinterpret_cast<char *>(&(tObj.mId)), sizeof(int));
-                
-                 if (samplesIn.eof())
-                        break;
-
-                Intersectable *sourceObj = NULL;
-                Intersectable *termObj = NULL;
-                
-                if (sObj.mId >= 0)
-                {
-                        ObjectContainer::iterator oit =
-                                lower_bound(objects.begin(), objects.end(), &sObj, ilt);
-                        sourceObj = *oit;
-                }
-                
-                if (tObj.mId >= 0)
-                {
-                        ObjectContainer::iterator oit =
-                                lower_bound(objects.begin(), objects.end(), &tObj, ilt);
-                        termObj = *oit;
-                }
-
-                samples.push_back(new VssRay(origin, termination, sourceObj, termObj));
-        }
-
-#endif
-        samplesIn.close();
-*/
-        return true;
-}
-
+                Plane3 plane = interior->GetPlane();
+                stream << "<Interior plane=\"" << plane.mNormal.x << " " 
+                           << plane.mNormal.y << " " << plane.mNormal.z << " " 
+                           << plane.mD << "\">" << endl;
+
+                ExportNode(interior->GetBack(), stream);
+                ExportNode(interior->GetFront(), stream);
+
+                stream << "</Interior>" << endl;
+        }
+}
 
 /************************************************************************/