source: GTP/trunk/Lib/Vis/Preprocessing/src/VspOspTree.cpp @ 1006

#include <stack>
#include <time.h>
#include <iomanip>

#include "Plane3.h"
#include "VspOspTree.h"
#include "Mesh.h"
#include "common.h"
#include "ViewCell.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include "Exporter.h"
#include "Plane3.h"
#include "ViewCellBsp.h"
#include "ViewCellsManager.h"
#include "Beam.h"

namespace GtpVisibilityPreprocessor {

#define USE_FIXEDPOINT_T 0


//-- static members

int VspOspTree::sFrontId = 0;
int VspOspTree::sBackId = 0;
int VspOspTree::sFrontAndBackId = 0;



// pvs penalty can be different from pvs size
inline static float EvalPvsPenalty(const int pvs, 
                                                                   const int lower,
                                                                   const int upper)
{
        // clamp to minmax values
        if (pvs < lower)
                return (float)lower;
        if (pvs > upper)
                return (float)upper;

        return (float)pvs;
}




/******************************************************************************/
/*                       class VspOspTree implementation                      */
/******************************************************************************/


VspOspTree::VspOspTree():
mRoot(NULL),
mOutOfBoundsCell(NULL),
mStoreRays(false),
mUseRandomAxis(false),
mTimeStamp(1)
{
        bool randomize = false;
        Environment::GetSingleton()->GetBoolValue("VspOspTree.Construction.randomize", randomize);
        if (randomize)
                Randomize(); // initialise random generator for heuristics

        //-- termination criteria for autopartition
        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.maxDepth", mTermMaxDepth);
        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.minPvs", mTermMinPvs);
        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.minRays", mTermMinRays);
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.minProbability", mTermMinProbability);
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxRayContribution", mTermMaxRayContribution);
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxCostRatio", mTermMaxCostRatio);
        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.missTolerance", mTermMissTolerance);
        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.maxViewCells", mMaxViewCells);

        //-- max cost ratio for early tree termination
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxCostRatio", mTermMaxCostRatio);

        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);

        // HACK//mTermMinPolygons = 25;

        //-- factors for bsp tree split plane heuristics
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.ct_div_ci", mCtDivCi);

        //-- partition criteria
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Construction.epsilon", mEpsilon);

        // if only the driving axis is used for axis aligned split
        Environment::GetSingleton()->GetBoolValue("VspOspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
        
        //Environment::GetSingleton()->GetFloatValue("VspOspTree.maxTotalMemory", mMaxTotalMemory);
        Environment::GetSingleton()->GetFloatValue("VspOspTree.maxStaticMemory", mMaxMemory);

        Environment::GetSingleton()->GetBoolValue("VspOspTree.useCostHeuristics", mUseCostHeuristics);
        Environment::GetSingleton()->GetBoolValue("VspOspTree.simulateOctree", mCirculatingAxis);
        Environment::GetSingleton()->GetBoolValue("VspOspTree.useRandomAxis", mUseRandomAxis);
        Environment::GetSingleton()->GetIntValue("VspOspTree.nodePriorityQueueType", mNodePriorityQueueType);
        
        char subdivisionStatsLog[100];
        Environment::GetSingleton()->GetStringValue("VspOspTree.subdivisionStats", subdivisionStatsLog);
        mSubdivisionStats.open(subdivisionStatsLog);

        Environment::GetSingleton()->GetFloatValue("VspOspTree.Construction.minBand", mMinBand);
        Environment::GetSingleton()->GetFloatValue("VspOspTree.Construction.maxBand", mMaxBand);
        

        //-- debug output

        Debug << "******* VSP BSP options ******** " << endl;
    Debug << "max depth: " << mTermMaxDepth << endl;
        Debug << "min PVS: " << mTermMinPvs << endl;
        Debug << "min probabiliy: " << mTermMinProbability << endl;
        Debug << "min rays: " << mTermMinRays << endl;
        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
        Debug << "miss tolerance: " << mTermMissTolerance << endl;
        Debug << "max view cells: " << mMaxViewCells << endl;
        Debug << "randomize: " << randomize << endl;

        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
        Debug << "max memory: " << mMaxMemory << endl;
        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
        Debug << "use random axis: " << mUseRandomAxis << endl;
        Debug << "priority queue type: " << mNodePriorityQueueType << endl;
        Debug << "circulating axis: " << mCirculatingAxis << endl;
        Debug << "minband: " << mMinBand << endl;
        Debug << "maxband: " << mMaxBand << endl;
        

        mSplitCandidates = new vector<SortableEntry>;

        Debug << endl;
}


BspViewCell *VspOspTree::GetOutOfBoundsCell()
{
        return mOutOfBoundsCell;
}


BspViewCell *VspOspTree::GetOrCreateOutOfBoundsCell()
{
        if (!mOutOfBoundsCell)
        {
                mOutOfBoundsCell = new BspViewCell();
                mOutOfBoundsCell->SetId(-1);
                mOutOfBoundsCell->SetValid(false);
        }

        return mOutOfBoundsCell;
}


const BspTreeStatistics &VspOspTree::GetStatistics() const
{
        return mBspStats;
}


VspOspTree::~VspOspTree()
{
        DEL_PTR(mRoot);
        DEL_PTR(mSplitCandidates);
}

void VspOspTree::Construct(const VssRayContainer &sampleRays,
                                                   AxisAlignedBox3 *forcedBoundingBox)
{
        mBspStats.nodes = 1;
        mBox.Initialize();      // initialise BSP tree bounding box

        if (forcedBoundingBox)
                mBox = *forcedBoundingBox;
        
        PolygonContainer polys;
        RayInfoContainer *rays = new RayInfoContainer();

        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();

        long startTime = GetTime();

        cout << "Extracting polygons from rays ... ";

        Intersectable::NewMail();

        int numObj = 0;

        //-- store rays
        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
        {
                VssRay *ray = *rit;

                float minT, maxT;

                static Ray hray;
                hray.Init(*ray);

                // TODO: not very efficient to implictly cast between rays types
                if (mBox.GetRaySegment(hray, minT, maxT))
                {
                        float len = ray->Length();

                        if (!len)
                                len = Limits::Small;

                        rays->push_back(RayInfo(ray, minT / len, maxT / len));
                }
        }

        mTermMinProbability *= mBox.GetVolume();

        mGlobalCostMisses = 0;

        cout << "finished" << endl;

        // use split cost priority queue
        Construct(rays);
}


// TODO: return memory usage in MB
float VspOspTree::GetMemUsage() const
{
        return (float)
                 (sizeof(VspOspTree) + 
                  mBspStats.Leaves() * sizeof(BspLeaf) + 
                  mCreatedViewCells * sizeof(BspViewCell) +
                  mBspStats.pvs * sizeof(ObjectPvsData) +
                  mBspStats.Interior() * sizeof(BspInterior) +
                  mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
}


void VspOspTree::Construct(RayInfoContainer *rays)
{
        VspOspSplitQueue tQueue;

        mRoot = new BspLeaf();

        const float prop = mBox.GetVolume();

        VspOspTraversalData tData(mRoot,
                                                          0,
                                                          rays,
                              ComputePvsSize(*rays),
                                                          prop,
                                                          mBox);


        // compute first split candidate
        VspOspSplitCandidate splitCandidate;
        EvalSplitCandidate(tData, splitCandidate);

        tQueue.push(splitCandidate);

        mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
        mTotalPvsSize = tData.mPvs;
        
        mSubdivisionStats 
                        << "#ViewCells\n1\n" <<  endl
                        << "#RenderCostDecrease\n0\n" << endl 
                        << "#SplitCandidateCost\n0\n" << endl
                        << "#TotalRenderCost\n" << mTotalCost << endl
                        << "#AvgRenderCost\n" << mTotalPvsSize << endl;

        Debug << "total cost: " << mTotalCost << endl;
        
        
        mBspStats.Start();
        cout << "Constructing vsp bsp tree ... \n";

        long startTime = GetTime();     
        int nLeaves = 500;
        int nViewCells = 500;

        // used for intermediate time measurements and progress
        long interTime = GetTime();     

        mOutOfMemory = false;

        mCreatedViewCells = 0;
        
        while (!tQueue.empty())
        {
                splitCandidate = tQueue.top();
            tQueue.pop();               

                // cost ratio of cost decrease / totalCost
                float costRatio = splitCandidate.GetCost() / mTotalCost;

                //Debug << "cost ratio: " << costRatio << endl;

                if (costRatio < mTermMinGlobalCostRatio)
                        ++ mGlobalCostMisses;
                
                if (0 && !mOutOfMemory)
                {
                        float mem = GetMemUsage();

                        if (mem > mMaxMemory)
                        {
                                mOutOfMemory = true;
                                Debug << "memory limit reached: " << mem << endl;
                        }
                }

                // subdivide leaf node
                BspNode *r = Subdivide(tQueue, splitCandidate);

                if (r == mRoot)
                        Debug << "VSP BSP tree construction time spent at root: "
                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;

                if (mBspStats.Leaves() >= nLeaves)
                {
                        nLeaves += 500;

                        cout << "leaves=" << mBspStats.Leaves() << endl;
                        Debug << "needed "
                                  << TimeDiff(interTime, GetTime())*1e-3 
                                  << " secs to create 500 view cells" << endl;
                        interTime = GetTime();
                }

                if (mCreatedViewCells == nViewCells)
                {
                        nViewCells += 500;

                        cout << "generated " << mCreatedViewCells << " viewcells" << endl;
                }
        }

        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
        cout << "finished\n";

        mBspStats.Stop();
}


bool VspOspTree::LocalTerminationCriteriaMet(const VspOspTraversalData &data) const
{
        return
                (((int)data.mRays->size() <= mTermMinRays) ||
                 (data.mPvs <= mTermMinPvs)   ||
                 (data.mProbability <= mTermMinProbability) ||
                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
                 (data.mDepth >= mTermMaxDepth));
}


bool VspOspTree::GlobalTerminationCriteriaMet(const VspOspTraversalData &data) const
{
        return
                (mOutOfMemory 
                || (mBspStats.Leaves() >= mMaxViewCells) 
                || (mGlobalCostMisses >= mTermGlobalCostMissTolerance)
                 );
}


// subdivide using a split plane queue
BspNode *VspOspTree::Subdivide(VspOspSplitQueue &tQueue,
                                                           VspOspSplitCandidate &splitCandidate)
{
        VspOspTraversalData &tData = splitCandidate.mParentData;

        BspNode *newNode = tData.mNode;

        if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
        {       
                VspOspTraversalData tFrontData;
                VspOspTraversalData tBackData;

                //-- continue subdivision
                
                // create new interior node and two leaf node
                //TODO
                const Plane3 splitPlane;// = splitCandidate.mSplitPlane;
                                
                newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData);
        
                const int maxCostMisses = splitCandidate.mMaxCostMisses;


                // how often was max cost ratio missed in this branch?
                tFrontData.mMaxCostMisses = maxCostMisses;
                tBackData.mMaxCostMisses = maxCostMisses;
                        
                
                if (1)
                {
                        float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
                        float cBack = (float)tBackData.mPvs * tBackData.mProbability;
                        float cData = (float)tData.mPvs * tData.mProbability;

                        
                        float costDecr = 
                                (cFront + cBack - cData) / mBox.GetVolume();

                        mTotalCost += costDecr;
                        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;

                        mSubdivisionStats 
                                        << "#ViewCells\n" << mBspStats.Leaves() << endl
                                        << "#RenderCostDecrease\n" << -costDecr << endl
                                        << "#SplitCandidateCost\n" << splitCandidate.GetCost() << endl
                                        << "#TotalRenderCost\n" << mTotalCost << endl
                                        << "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mBspStats.Leaves() << endl;
                }

        
                //-- push the new split candidates on the stack
                VspOspSplitCandidate frontCandidate;
                VspOspSplitCandidate backCandidate;

                EvalSplitCandidate(tFrontData, frontCandidate);
                EvalSplitCandidate(tBackData, backCandidate);
        
                tQueue.push(frontCandidate);
                tQueue.push(backCandidate);
        
                // delete old leaf node
                DEL_PTR(tData.mNode);
        }


        //-- terminate traversal and create new view cell
        if (newNode->IsLeaf())
        {
                BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);
                BspViewCell *viewCell = new BspViewCell();

        leaf->SetViewCell(viewCell);
                
                //-- update pvs
                int conSamp = 0;
                float sampCon = 0.0f;
                AddToPvs(leaf, *tData.mRays, sampCon, conSamp);

                // update scalar pvs size value
                mViewCellsManager->SetScalarPvsSize(viewCell, viewCell->GetPvs().GetSize());

                mBspStats.contributingSamples += conSamp;
                mBspStats.sampleContributions +=(int) sampCon;

                //-- store additional info
                if (mStoreRays)
                {
                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
                        for (it = tData.mRays->begin(); it != it_end; ++ it)
                        {
                                (*it).mRay->Ref();                      
                                leaf->mVssRays.push_back((*it).mRay);
                        }
                }

                // should I check here?
                viewCell->mLeaf = leaf;

                viewCell->SetVolume(tData.mProbability);
        leaf->mProbability = tData.mProbability;

                EvaluateLeafStats(tData);               
        }

        //-- cleanup
        tData.Clear();

        return newNode;
}


void VspOspTree::EvalSplitCandidate(VspOspTraversalData &tData,
                                                                        VspOspSplitCandidate &splitData)
{
        VspOspTraversalData frontData;
        VspOspTraversalData backData;

        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);

        // compute locally best split plane
        // TODO
    const bool success =0;/// SelectPlane(splitData.mSplitPlane, tData, 
//                                                                       frontData, backData);
        //TODO
        // compute global decrease in render cost
        splitData.mRenderCost = 0;//EvalRenderCostDecrease(splitData.mSplitPlane, tData);
        splitData.mParentData = tData;
        splitData.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
}


BspInterior *VspOspTree::SubdivideNode(const Plane3 &splitPlane,
                                                                           VspOspTraversalData &tData,
                                                                           VspOspTraversalData &frontData,
                                                                           VspOspTraversalData &backData)
{
        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
        
        //-- the front and back traversal data is filled with the new values
        frontData.mDepth = tData.mDepth + 1;
        frontData.mRays = new RayInfoContainer();
        
        backData.mDepth = tData.mDepth + 1;
        backData.mRays = new RayInfoContainer();
        

        //-- subdivide rays
        SplitRays(splitPlane,
                          *tData.mRays,
                          *frontData.mRays,
                          *backData.mRays);


        // compute pvs
        frontData.mPvs = ComputePvsSize(*frontData.mRays);
        backData.mPvs = ComputePvsSize(*backData.mRays);

        // split front and back node geometry and compute area
        
        // if geometry was not already computed
        // TODO

#if TODO
        tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition, 
                                                         frontData.mBoundingBox, backData.mBoundingBox);

#endif

        frontData.mProbability = frontData.mBoundingBox.GetVolume();
        backData.mProbability = tData.mProbability - frontData.mProbability;

        
    ///////////////////////////////////////////
        // subdivide further

        // store maximal and minimal depth
        if (tData.mDepth > mBspStats.maxDepth)
        {
                Debug << "max depth increases to " << tData.mDepth << " at " << mBspStats.Leaves() << " leaves" << endl;
                mBspStats.maxDepth = tData.mDepth;
        }

        mBspStats.nodes += 2;

    
        BspInterior *interior = new BspInterior(splitPlane);

#ifdef _DEBUG
        Debug << interior << endl;
#endif


        //-- create front and back leaf

        BspInterior *parent = leaf->GetParent();

        // replace a link from node's parent
        if (parent)
        {
                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();

        interior->mTimeStamp = mTimeStamp ++;
        
        return interior;
}


void VspOspTree::AddToPvs(BspLeaf *leaf,
                                                  const RayInfoContainer &rays,
                                                  float &sampleContributions,
                                                  int &contributingSamples)
{
        sampleContributions = 0;
        contributingSamples = 0;
  
        RayInfoContainer::const_iterator it, it_end = rays.end();
  
        ViewCellLeaf *vc = leaf->GetViewCell();
  
        // add contributions from samples to the PVS
        for (it = rays.begin(); it != it_end; ++ it)
        {
                float sc = 0.0f;
                VssRay *ray = (*it).mRay;

                bool madeContrib = false;
                float contribution;

                if (ray->mTerminationObject) 
                {
                        if (vc->AddPvsSample(ray->mTerminationObject, ray->mPdf, contribution))
                                madeContrib = true;
                        sc += contribution;
                }
          
                if (ray->mOriginObject) 
                {
                        if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
                                madeContrib = true;
                        sc += contribution;
                }
          
          sampleContributions += sc;

          if (madeContrib)
                  ++ contributingSamples;
                
          //leaf->mVssRays.push_back(ray);
        }
}


void VspOspTree::SortSplitCandidates(const RayInfoContainer &rays, 
                                                                         const int axis, 
                                                                         float minBand, 
                                                                         float maxBand)
{
        mSplitCandidates->clear();

        int requestedSize = 2 * (int)(rays.size());
        // creates a sorted split candidates array
        if (mSplitCandidates->capacity() > 500000 &&
                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
        {
        delete mSplitCandidates;
                mSplitCandidates = new vector<SortableEntry>;
        }

        mSplitCandidates->reserve(requestedSize);

        float pos;

        // float values => don't compare with exact values
        if (0)
        {
                minBand += Limits::Small;
                maxBand -= Limits::Small;
        }

        // insert all queries
        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
        {
                const bool positive = (*ri).mRay->HasPosDir(axis);
                                
                pos = (*ri).ExtrapOrigin(axis);
                // clamp to min / max band
                if (0) ClipValue(pos, minBand, maxBand);
                
                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax, 
                                                                        pos, (*ri).mRay));

                pos = (*ri).ExtrapTermination(axis);
                // clamp to min / max band
                if (0) ClipValue(pos, minBand, maxBand);

                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin, 
                                                                        pos, (*ri).mRay));
        }

        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
}


float VspOspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
                                                                                  const AxisAlignedBox3 &box,
                                                                                  const int pvsSize,
                                                                                  const int axis,
                                          float &position)
{
        const float minBox = box.Min(axis);
        const float maxBox = box.Max(axis);

        const float sizeBox = maxBox - minBox;

        const float minBand = minBox + mMinBand * sizeBox;
        const float maxBand = minBox + mMaxBand * sizeBox;

        SortSplitCandidates(rays, axis, minBand, maxBand);

        // go through the lists, count the number of objects left and right
        // and evaluate the following cost funcion:
        // C = ct_div_ci  + (ql*rl + qr*rr)/queries

        int pvsl = 0;
        int pvsr = pvsSize;

        int pvsBack = pvsl;
        int pvsFront = pvsr;

        float sum = (float)pvsSize * sizeBox;
        float minSum = 1e20f;

        
        // if no border can be found, take mid split
        position = minBox + 0.5f * sizeBox;
        
        // the relative cost ratio
        float ratio = /*Limits::Infinity;*/99999999.0f;
        bool splitPlaneFound = false;

        Intersectable::NewMail();

        RayInfoContainer::const_iterator ri, ri_end = rays.end();

        // set all object as belonging to the front pvs
        for(ri = rays.begin(); ri != ri_end; ++ ri)
        {
                Intersectable *oObject = (*ri).mRay->mOriginObject;
                Intersectable *tObject = (*ri).mRay->mTerminationObject;

                if (oObject)
                {
                        if (!oObject->Mailed())
                        {
                                oObject->Mail();
                                oObject->mCounter = 1;
                        }
                        else
                        {
                                ++ oObject->mCounter;
                        }
                }
                if (tObject)
                {
                        if (!tObject->Mailed())
                        {
                                tObject->Mail();
                                tObject->mCounter = 1;
                        }
                        else
                        {
                                ++ tObject->mCounter;
                        }
                }
        }

        Intersectable::NewMail();

        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();

        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
        {
                VssRay *ray;
                ray = (*ci).ray;
                
                Intersectable *oObject = ray->mOriginObject;
                Intersectable *tObject = ray->mTerminationObject;
                

                switch ((*ci).type)
                {
                        case SortableEntry::ERayMin:
                                {
                                        if (oObject && !oObject->Mailed())
                                        {
                                                oObject->Mail();
                                                ++ pvsl;
                                        }
                                        if (tObject && !tObject->Mailed())
                                        {
                                                tObject->Mail();
                                                ++ pvsl;
                                        }
                                        break;
                                }
                        case SortableEntry::ERayMax:
                                {
                                        if (oObject)
                                        {
                                                if (-- oObject->mCounter == 0)
                                                        -- pvsr;
                                        }

                                        if (tObject)
                                        {
                                                if (-- tObject->mCounter == 0)
                                                        -- pvsr;
                                        }

                                        break;
                                }
                }

                
                
                // Note: sufficient to compare size of bounding boxes of front and back side?
                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
                {
                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);

                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << endl;
                        //Debug << "cost= " << sum << endl;

                        if (sum < minSum)
                        {
                                splitPlaneFound = true;

                                minSum = sum;
                                position = (*ci).value;
                                
                                pvsBack = pvsl;
                                pvsFront = pvsr;
                        }
                }
        }
        
        
        // -- compute cost
        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();

        const float pOverall = sizeBox;

        const float pBack = position - minBox;
        const float pFront = maxBox - position;
        
        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
        
        const float oldRenderCost = penaltyOld * pOverall;
        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

        if (splitPlaneFound)
        {
                ratio = newRenderCost / (oldRenderCost + Limits::Small);
        }
        //if (axis != 1)
        //Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
         //    <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;

        return ratio;
}


float VspOspTree::SelectPlane(const VspOspTraversalData &tData,
                                                          AxisAlignedPlane &plane,
                                                          float &pFront,
                                                          float &pBack)
{
        float nPosition[3];
        float nCostRatio[3];
        float nProbFront[3];
        float nProbBack[3];

        // create bounding box of node geometry
        AxisAlignedBox3 box = tData.mBoundingBox;
                
        int sAxis = 0;
        int bestAxis = -1;

#if 0
        // maximum cost ratio for axis to be valid:
        // if exceeded, spatial mid split is used instead
        const maxCostRatioForHeur = 0.99f;
#endif

        // if we use some kind of specialised fixed axis
    const bool useSpecialAxis = 
                mOnlyDrivingAxis || mUseRandomAxis || mCirculatingAxis;

        if (mUseRandomAxis)
                sAxis = Random(3);
        else if (mCirculatingAxis)
                sAxis = (tData.mAxis + 1) % 3;
        else if (mOnlyDrivingAxis)
                sAxis = box.Size().DrivingAxis();


        for (int axis = 0; axis < 3 ; ++ axis)
        {
                if (!useSpecialAxis || (axis == sAxis))
                {
                        //-- place split plane using heuristics

                        if (mUseCostHeuristics)
                        {
                                nCostRatio[axis] =
                                        BestCostRatioHeuristics(*tData.mRays,
                                                                                    box,
                                                                                        tData.mPvs,
                                                                                        axis,
                                                                                        nPosition[axis]);                       
                        }
                        else //-- split plane position is spatial median
                        {
                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;

                                nCostRatio[axis] = EvalSplitCost(tData,
                                                                                                 box,
                                                                                                 axis,
                                                                                                 nPosition[axis],
                                                                                                 nProbFront[axis], 
                                                                                                 nProbBack[axis]);
                        }
                                                
                        if (bestAxis == -1)
                        {
                                bestAxis = axis;
                        }
                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
                        {
                                bestAxis = axis;
                        } 
                }
        }

        //-- assign values
        plane.mAxis = bestAxis;
        pFront = nProbFront[bestAxis];
        pBack = nProbBack[bestAxis];

        //-- split plane position
    plane.mPosition = nPosition[bestAxis];

        return nCostRatio[bestAxis];
}


inline void VspOspTree::GenerateUniqueIdsForPvs()
{
        Intersectable::NewMail(); sBackId = Intersectable::sMailId;
        Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
        Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
}


float VspOspTree::EvalRenderCostDecrease(const Plane3 &candidatePlane,
                                                                                 const VspOspTraversalData &data) const
{
        float pvsFront = 0;
        float pvsBack = 0;
        float totalPvs = 0;

        // probability that view point lies in back / front node
        float pOverall = data.mProbability;
        float pFront = 0;
        float pBack = 0;


        // create unique ids for pvs heuristics
        GenerateUniqueIdsForPvs();
        
        for (int i = 0; i < data.mRays->size(); ++ i)
        {
                RayInfo rayInf = (*data.mRays)[i];

                float t;
                VssRay *ray = rayInf.mRay;
                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);

                // find front and back pvs for origing and termination object
                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
        }


        BspNodeGeometry geomFront;
        BspNodeGeometry geomBack;
        
        pFront = geomFront.GetVolume();
        pBack = pOverall - pFront;

                

        // -- pvs rendering heuristics
        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();

        // only render cost heuristics or combined with standard deviation
        const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
                        
        const float oldRenderCost = pOverall * penaltyOld;
        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;

        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
        const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBox.GetVolume();
        

#if 1
        // take render cost of node into account to avoid being stuck in a local minimum
        const float normalizedOldRenderCost = oldRenderCost / mBox.GetVolume();
        return 0.99f * renderCostDecrease + 0.01f * normalizedOldRenderCost;
#else
        return renderCostDecrease;
#endif
}



float VspOspTree::EvalSplitCost(const VspOspTraversalData &data,
                                                                const AxisAlignedBox3 &box,
                                                                const int axis,
                                                                const float &position,                                                                            
                                                                float &pFront,
                                                                float &pBack) const
{
        float pvsTotal = 0;
        float pvsFront = 0;
        float pvsBack = 0;
        
        // create unique ids for pvs heuristics
        GenerateUniqueIdsForPvs();

        const int pvsSize = data.mPvs;

        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();

        // this is the main ray classification loop!
        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
        {
                // determine the side of this ray with respect to the plane
                float t;
                const int side = (*rit).ComputeRayIntersection(axis, position, t);
        
                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
        }

        //-- pvs heuristics
        float pOverall;

        //-- compute heurstics
        //-- we take simplified computation for mid split
                
        pOverall = data.mProbability;
        pBack = pFront = pOverall * 0.5f;
        
        
#ifdef _DEBUG
        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
        Debug << pFront << " " << pBack << " " << pOverall << endl;
#endif

        
        const float newCost = pvsBack * pBack + pvsFront * pFront;
        const float oldCost = (float)pvsSize * pOverall + Limits::Small;

        return  (mCtDivCi + newCost) / oldCost;
}


void VspOspTree::AddObjToPvs(Intersectable *obj,
                                                         const int cf,
                                                         float &frontPvs,
                                                         float &backPvs,
                                                         float &totalPvs) const
{
        if (!obj)
                return;

        const float renderCost = mViewCellsManager->EvalRenderCost(obj);

        // new object
        if ((obj->mMailbox != sFrontId) &&
                (obj->mMailbox != sBackId) &&
                (obj->mMailbox != sFrontAndBackId))
        {
                totalPvs += renderCost;
        }

        // 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 += renderCost;
                
                        if (obj->mMailbox == sBackId)
                                obj->mMailbox = sFrontAndBackId;
                        else
                                obj->mMailbox = sFrontId;
                }
        }

        if (cf <= 0)
        {
                if ((obj->mMailbox != sBackId) &&
                        (obj->mMailbox != sFrontAndBackId))
                {
                        backPvs += renderCost;
                
                        if (obj->mMailbox == sFrontId)
                                obj->mMailbox = sFrontAndBackId;
                        else
                                obj->mMailbox = sBackId;
                }
        }
}


void VspOspTree::CollectLeaves(vector<BspLeaf *> &leaves, 
                                                           const bool onlyUnmailed,
                                                           const int maxPvsSize) const
{
        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);

        while (!nodeStack.empty())
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
                
                if (node->IsLeaf())
                {
                        // test if this leaf is in valid view space
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
                        if (leaf->TreeValid() && 
                                (!onlyUnmailed || !leaf->Mailed()) &&
                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
                        {
                                leaves.push_back(leaf);
                        }
                }
                else
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

                        nodeStack.push(interior->GetBack());
                        nodeStack.push(interior->GetFront());
                }
        }
}


AxisAlignedBox3 VspOspTree::GetBoundingBox() const
{
        return mBox;
}


BspNode *VspOspTree::GetRoot() const
{
        return mRoot;
}


void VspOspTree::EvaluateLeafStats(const VspOspTraversalData &data)
{
        // the node became a leaf -> evaluate stats for leafs
        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);


        if (data.mPvs > mBspStats.maxPvs)
        {
                mBspStats.maxPvs = data.mPvs;
        }

        mBspStats.pvs += data.mPvs;

        if (data.mDepth < mBspStats.minDepth)
        {
                mBspStats.minDepth = data.mDepth;
        }
        
        if (data.mDepth >= mTermMaxDepth)
        {
        ++ mBspStats.maxDepthNodes;
                //Debug << "new max depth: " << mBspStats.maxDepthNodes << endl;
        }

        // accumulate rays to compute rays /  leaf
        mBspStats.accumRays += (int)data.mRays->size();

        if (data.mPvs < mTermMinPvs)
                ++ mBspStats.minPvsNodes;

        if ((int)data.mRays->size() < mTermMinRays)
                ++ mBspStats.minRaysNodes;

        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
                ++ mBspStats.maxRayContribNodes;

        if (data.mProbability <= mTermMinProbability)
                ++ mBspStats.minProbabilityNodes;
        
        // accumulate depth to compute average depth
        mBspStats.accumDepth += data.mDepth;

        ++ mCreatedViewCells;

#ifdef _DEBUG
        Debug << "BSP stats: "
                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
          //              << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
                  << "#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 VspOspTree::CastRay(Ray &ray)
{
        int hits = 0;

        stack<BspRayTraversalData> tQueue;

        float maxt, mint;

        if (!mBox.GetRaySegment(ray, mint, maxt))
                return 0;

        Intersectable::NewMail();
        ViewCell::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
                        tQueue.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::VspOspIntersection(maxt, leaf));
                                leaf->GetViewCell()->Mail();
                                ++ hits;
                        }

                        //-- fetch the next far child from the stack
                        if (tQueue.empty())
                                break;

                        entp = extp;
                        mint = maxt; // NOTE: need this?

                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
                                break;

                        BspRayTraversalData &s = tQueue.top();

                        node = s.mNode;
                        extp = s.mExitPoint;
                        maxt = s.mMaxT;

                        tQueue.pop();
                }
        }

        return hits;
}


void VspOspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
{
        ViewCell::NewMail();
        
        CollectViewCells(mRoot, onlyValid, viewCells, true);
}


void VspOspTree::CollapseViewCells()
{
// TODO
#if HAS_TO_BE_REDONE
        stack<BspNode *> nodeStack;

        if (!mRoot)
                return;

        nodeStack.push(mRoot);
        
        while (!nodeStack.empty()) 
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
                
                if (node->IsLeaf())
        {
                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();

                        if (!viewCell->GetValid())
                        {
                                BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
        
                                ViewCellContainer leaves;
                                mViewCellsTree->CollectLeaves(viewCell, leaves);

                                ViewCellContainer::const_iterator it, it_end = leaves.end();

                                for (it = leaves.begin(); it != it_end; ++ it)
                                {
                                        BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
                                        ++ mBspStats.invalidLeaves;
                                }

                                // add to unbounded view cell
                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
                                DEL_PTR(viewCell);
                        }
                }
                else
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
                
                        nodeStack.push(interior->GetFront());
                        nodeStack.push(interior->GetBack());
                }
        }

        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
#endif
}


void VspOspTree::CollectRays(VssRayContainer &rays)
{
        vector<BspLeaf *> leaves;

        vector<BspLeaf *>::const_iterator lit, lit_end = leaves.end();

        for (lit = leaves.begin(); lit != lit_end; ++ lit)
        {
                BspLeaf *leaf = *lit;
                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();

                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
                        rays.push_back(*rit);
        }
}


void VspOspTree::ValidateTree()
{
        stack<BspNode *> nodeStack;

        if (!mRoot)
                return;

        nodeStack.push(mRoot);
        
        mBspStats.invalidLeaves = 0;
        while (!nodeStack.empty()) 
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
                
                if (node->IsLeaf())
                {
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);

                        if (!leaf->GetViewCell()->GetValid())
                                ++ mBspStats.invalidLeaves;

                        // validity flags don't match => repair
                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
                        {
                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
                                PropagateUpValidity(leaf);
                        }
                }
                else
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
                
                        nodeStack.push(interior->GetFront());
                        nodeStack.push(interior->GetBack());
                }
        }

        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
}



void VspOspTree::CollectViewCells(BspNode *root,
                                                                  bool onlyValid,
                                                                  ViewCellContainer &viewCells,
                                                                  bool onlyUnmailed) const
{
        stack<BspNode *> nodeStack;

        if (!root)
                return;

        nodeStack.push(root);
        
        while (!nodeStack.empty()) 
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
                
                if (node->IsLeaf())
                {
                        if (!onlyValid || node->TreeValid())
                        {
                                ViewCellLeaf *leafVc = dynamic_cast<BspLeaf *>(node)->GetViewCell();

                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
                                                
                                if (!onlyUnmailed || !viewCell->Mailed()) 
                                {
                                        viewCell->Mail();
                                        viewCells.push_back(viewCell);
                                }
                        }
                }
                else
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
                
                        nodeStack.push(interior->GetFront());
                        nodeStack.push(interior->GetBack());
                }
        }

}


int VspOspTree::SplitRays(const Plane3 &plane,
                                                  RayInfoContainer &rays,
                                                  RayInfoContainer &frontRays,
                                                  RayInfoContainer &backRays) const
{
        int splits = 0;

        RayInfoContainer::const_iterator it, it_end = rays.end();

        for (it = rays.begin(); it != it_end; ++ it)
        {
                RayInfo bRay = *it;
                
                VssRay *ray = bRay.mRay;
                float t;

                // get classification and receive new t
                const int cf = bRay.ComputeRayIntersection(plane, t);

                switch (cf)
                {
                case -1:
                        backRays.push_back(bRay);
                        break;
                case 1:
                        frontRays.push_back(bRay);
                        break;
                case 0:
                        {
                                //-- split ray
                                //-- test if start point behind or in front of plane
                                const int side = plane.Side(bRay.ExtrapOrigin());

                                if (side <= 0)
                                {
                                        backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
                                        frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
                                }
                                else
                                {
                                        frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
                                        backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
                                }
                        }
                        break;
                default:
                        Debug << "Should not come here" << endl;
                        break;
                }
        }

        return splits;
}


int VspOspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
                                                          const bool onlyUnmailed) const
{
        // TODO
        return 0;
}


BspLeaf *VspOspTree::GetRandomLeaf(const Plane3 &halfspace)
{
#if TODO
    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;
                        BspNodeGeometry geom;

                        // todo: not very efficient: constructs full cell everytime
                        ConstructGeometry(interior, geom);

                        const int cf = 
                                Polygon3::ClassifyPlane(geom.GetPolys(), 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);
                }
        }
#endif
        return NULL;
}


BspLeaf *VspOspTree::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 VspOspTree::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 VspOspTree::GetEpsilon() const
{
        return mEpsilon;
}



int VspOspTree::CastLineSegment(const Vector3 &origin,
                                                                const Vector3 &termination,
                                                                ViewCellContainer &viewcells)
{
        int hits = 0;
        stack<BspRayTraversalData> tStack;

        float mint = 0.0f, maxt = 1.0f;

        Intersectable::NewMail();
        ViewCell::NewMail();

        Vector3 entp = origin;
        Vector3 extp = termination;

        BspNode *node = mRoot;
        BspNode *farChild = NULL;

        float t;
        const float thresh = 1e-6f; // matt: change this to adjustable value
        
        while (1)
        {
                if (!node->IsLeaf())
                {
                        BspInterior *in = dynamic_cast<BspInterior *>(node);

                        Plane3 splitPlane = in->GetPlane();
                        
                        const int entSide = splitPlane.Side(entp, thresh);
                        const int extSide = splitPlane.Side(extp, thresh);

                        if (entSide < 0) 
                        {
                                node = in->GetBack();
                                
                                // plane does not split ray => no far child
                                if (extSide <= 0) 
                                        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 // one of the ray end points is on the plane
                        {       // NOTE: what to do if ray is coincident with plane?
                                if (extSide < 0)
                                        node = in->GetBack();
                                else //if (extSide > 0)
                                        node = in->GetFront();
                                //else break; // coincident => count no intersections

                                continue; // no far child
                        }

                        // push data for far child
                        tStack.push(BspRayTraversalData(farChild, extp));

                        // find intersection of ray segment with plane
                        extp = splitPlane.FindIntersection(origin, extp, &t);
                } 
                else
                {
                        // reached leaf => intersection with view cell
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
                        ViewCell *viewCell;
                        
                        if (0)
                                viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
                        else
                                viewCell = leaf->GetViewCell();

                        if (!viewCell->Mailed())
                        {
                                viewcells.push_back(viewCell);
                                viewCell->Mail();
                                ++ hits;
                        }

                        //-- fetch the next far child from the stack
                        if (tStack.empty())
                                break;

                        entp = extp;
                        
                        const BspRayTraversalData &s = tStack.top();

                        node = s.mNode;
                        extp = s.mExitPoint;

                        tStack.pop();
                }
        }

        return hits;
}



BspNode *VspOspTree::CollapseTree(BspNode *node, int &collapsed)
{
// TODO
#if HAS_TO_BE_REDONE
        if (node->IsLeaf())
                return node;

        BspInterior *interior = dynamic_cast<BspInterior *>(node);

        BspNode *front = CollapseTree(interior->GetFront(), collapsed);
        BspNode *back = CollapseTree(interior->GetBack(), collapsed);

        if (front->IsLeaf() && back->IsLeaf())
        {
                BspLeaf *frontLeaf = dynamic_cast<BspLeaf *>(front);
                BspLeaf *backLeaf = dynamic_cast<BspLeaf *>(back);

                //-- collapse tree
                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
                {
                        BspViewCell *vc = frontLeaf->GetViewCell();

                        BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);
                        leaf->SetTreeValid(frontLeaf->TreeValid());

                        // replace a link from node's parent
                        if (leaf->GetParent())
                                leaf->GetParent()->ReplaceChildLink(node, leaf);
                        else
                                mRoot = leaf;

                        ++ collapsed;
                        delete interior;

                        return leaf;
                }
        }
#endif
        return node;
}


int VspOspTree::CollapseTree()
{
        int collapsed = 0;
        //TODO
#if HAS_TO_BE_REDONE
        (void)CollapseTree(mRoot, collapsed);

        // revalidate leaves
        RepairViewCellsLeafLists();
#endif
        return collapsed;
}


void VspOspTree::RepairViewCellsLeafLists()
{
// TODO
#if HAS_TO_BE_REDONE
        // list not valid anymore => clear
        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);

        ViewCell::NewMail();

        while (!nodeStack.empty())
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();

                if (node->IsLeaf())
                {
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);

                        BspViewCell *viewCell = leaf->GetViewCell();

                        if (!viewCell->Mailed())
                        {
                                viewCell->mLeaves.clear();
                                viewCell->Mail();
                        }
        
                        viewCell->mLeaves.push_back(leaf);

                }
                else
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

                        nodeStack.push(interior->GetFront());
                        nodeStack.push(interior->GetBack());
                }
        }
// TODO
#endif
}



ViewCell *VspOspTree::GetViewCell(const Vector3 &point, const bool active)
{
        if (mRoot == NULL)
                return NULL;

        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);
  
        ViewCellLeaf *viewcell = NULL;
  
        while (!nodeStack.empty())  
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
        
                if (node->IsLeaf()) 
                {
                        viewcell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
                        break;
                } 
                else    
                {       
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
        
                        // random decision
                        if (interior->GetPlane().Side(point) < 0)
                                nodeStack.push(interior->GetBack());
                        else
                                nodeStack.push(interior->GetFront());
                }
        }
  
        if (active)
                return mViewCellsTree->GetActiveViewCell(viewcell);
        else
                return viewcell;
}


bool VspOspTree::ViewPointValid(const Vector3 &viewPoint) const
{
        BspNode *node = mRoot;

        while (1)
        {
                // early exit
                if (node->TreeValid())
                        return true;

                if (node->IsLeaf())
                        return false;
                        
                BspInterior *in = dynamic_cast<BspInterior *>(node);
                                        
                if (in->GetPlane().Side(viewPoint) <= 0) 
                {
                        node = in->GetBack();
                } 
                else
                {
                        node = in->GetFront();
                }
        }

        // should never come here
        return false;
}


void VspOspTree::PropagateUpValidity(BspNode *node)
{
        const bool isValid = node->TreeValid();

        // propagative up invalid flag until only invalid nodes exist over this node
        if (!isValid)
        {
                while (!node->IsRoot() && node->GetParent()->TreeValid())
                {
                        node = node->GetParent();
                        node->SetTreeValid(false);
                }
        }
        else
        {
                // propagative up valid flag until one of the subtrees is invalid
                while (!node->IsRoot() && !node->TreeValid())
                {
            node = node->GetParent();
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
                        
                        // the parent is valid iff both leaves are valid
                        node->SetTreeValid(interior->GetBack()->TreeValid() && 
                                                           interior->GetFront()->TreeValid());
                }
        }
}

#if ZIPPED_VIEWCELLS
bool VspOspTree::Export(ogzstream &stream)
#else
bool VspOspTree::Export(ofstream &stream)
#endif
{
        ExportNode(mRoot, stream);

        return true;
}


#if ZIPPED_VIEWCELLS
void VspOspTree::ExportNode(BspNode *node, ogzstream &stream)
#else
void VspOspTree::ExportNode(BspNode *node, ofstream &stream)
#endif
{
        if (node->IsLeaf())
        {
                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());

                int id = -1;
                if (viewCell != mOutOfBoundsCell)
                        id = viewCell->GetId();

                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
        }
        else
        {
                BspInterior *interior = dynamic_cast<BspInterior *>(node);
        
                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;
        }
}

}