source: GTP/trunk/Lib/Vis/Preprocessing/src/BvHierarchy.cpp @ 1251

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

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


namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0

static float debugVol = 0;
int BvhNode::sMailId = 0;
BvHierarchy *BvHierarchy::BvhSubdivisionCandidate::sBvHierarchy = NULL;


inline static bool ilt(Intersectable *obj1, Intersectable *obj2)
{
        return obj1->mId < obj2->mId;
}


/***************************************************************/
/*              class BvhNode implementation                   */
/***************************************************************/

BvhNode::BvhNode(): mParent(NULL)
{
}

BvhNode::BvhNode(const AxisAlignedBox3 &bbox):
mBoundingBox(bbox)
{
}


BvhNode::BvhNode(const AxisAlignedBox3 &bbox, BvhInterior *parent):
mBoundingBox(bbox), mParent(parent)
{
}


bool BvhNode::IsRoot() const
{
        return mParent == NULL;
}


BvhInterior *BvhNode::GetParent()
{
        return mParent;
}


void BvhNode::SetParent(BvhInterior *parent)
{
        mParent = parent;
}



/******************************************************************/
/*              class BvhInterior implementation                  */
/******************************************************************/


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox):
BvhNode(bbox)
{
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent):
BvhNode(bbox, parent)
{
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent, const int numObjects):
BvhNode(bbox, parent)
{
        mObjects.reserve(numObjects);
}


bool BvhLeaf::IsLeaf() const 
{ 
        return true; 
}


BvhLeaf::~BvhLeaf() 
{ 
}


/******************************************************************/
/*              class BvhInterior implementation                  */
/******************************************************************/


BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox):
BvhNode(bbox)
{
}


BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox, BvhInterior *parent):
BvhNode(bbox, parent)
{
}


void BvhInterior::ReplaceChildLink(BvhNode *oldChild, BvhNode *newChild) 
{
        if (mBack == oldChild)
                mBack = newChild;
        else
                mFront = newChild;
}


bool BvhInterior::IsLeaf() const 
{ 
        return false; 
}


BvhInterior::~BvhInterior() 
{ 
        DEL_PTR(mFront); 
        DEL_PTR(mBack);
}


void BvhInterior::SetupChildLinks(BvhNode *front, BvhNode *back) 
{
    mBack = back;
    mFront = front;
}



/*******************************************************************/
/*                  class BvHierarchy implementation               */
/*******************************************************************/


BvHierarchy::BvHierarchy():
mRoot(NULL),
mTimeStamp(1)
{
        ReadEnvironment();
        mSubdivisionCandidates = new vector<SortableEntry>;
}


BvHierarchy::~BvHierarchy()
{
        // delete kd intersectables
        BvhIntersectableMap::iterator it, it_end = mBvhIntersectables.end();

        for (it = mBvhIntersectables.begin(); it != mBvhIntersectables.end(); ++ it)
        {
                DEL_PTR((*it).second);
        }

        DEL_PTR(mSubdivisionCandidates);

        mSubdivisionStats.close();
}


void BvHierarchy::ReadEnvironment()
{
        bool randomize = false;
        Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
        if (randomize)
                Randomize(); // initialise random generator for heuristics

        //-- termination criteria for autopartition
        Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxDepth", mTermMaxDepth);
        Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxLeaves", mTermMaxLeaves);
        Environment::GetSingleton()->GetIntValue("OspTree.Termination.minObjects", mTermMinObjects);
        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minProbability", mTermMinVolume);
        
        Environment::GetSingleton()->GetIntValue("OspTree.Termination.missTolerance", mTermMissTolerance);
        
        //-- max cost ratio for early tree termination
        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.maxCostRatio", mTermMaxCostRatio);

        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minGlobalCostRatio",
                mTermMinGlobalCostRatio);
        

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

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

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


        char subdivisionStatsLog[100];
        Environment::GetSingleton()->GetStringValue("OspTree.subdivisionStats", subdivisionStatsLog);
        mSubdivisionStats.open(subdivisionStatsLog);

        Environment::GetSingleton()->GetFloatValue("OspTree.Construction.splitBorder", mSplitBorder);
        Environment::GetSingleton()->GetFloatValue("OspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
        

        //-- debug output

        Debug << "******* Bvh hierarchy options ******** " << endl;

    Debug << "max depth: " << mTermMaxDepth << endl;
        Debug << "min probabiliy: " << mTermMinVolume<< endl;
        Debug << "min objects: " << mTermMinObjects << endl;
        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
        Debug << "miss tolerance: " << mTermMissTolerance << endl;
        Debug << "max leaves: " << mTermMaxLeaves << 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 << "split borders: " << mSplitBorder << endl;
        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
        Debug << endl;
}


void AssociateObjectsWithLeaf(BvhLeaf *leaf)
{
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                (*oit)->mBvhLeaf = leaf;
        }
}


BvhInterior *BvHierarchy::SubdivideNode(const ObjectContainer &frontObjects,
                                                                                const ObjectContainer &backObjects,
                                                                                const BvhTraversalData &tData,
                                                                                BvhTraversalData &frontData,
                                                                                BvhTraversalData &backData)
{
        BvhLeaf *leaf = tData.mNode;
        
        // two new leaves
    mBvhStats.nodes += 2;

        // add the new nodes to the tree
        BvhInterior *node = new BvhInterior(tData.mBoundingBox, leaf->GetParent());


        //-- the front and back traversal data is filled with the new values

        frontData.mDepth = backData.mDepth = tData.mDepth + 1;
                
        // frontData.mRays = new RayInfoContainer();
        // backData.mRays = new RayInfoContainer();

        frontData.mBoundingBox = ComputeBoundingBox(frontObjects);
        backData.mBoundingBox = ComputeBoundingBox(backObjects);
        //RemoveParentViewCellsPvs(leaf, *tData.mRays);


        /////////////
        //-- create front and back leaf

        BvhLeaf *back = new BvhLeaf(backData.mBoundingBox, node, (int)backObjects.size());
        BvhLeaf *front = new BvhLeaf(frontData.mBoundingBox, node, (int)frontObjects.size());

        BvhInterior *parent = leaf->GetParent();


        // replace a link from node's parent
        if (parent)
        {
                parent->ReplaceChildLink(leaf, node);
                node->SetParent(parent);
        }
        else // new root
        {
                mRoot = node;
        }

        // and setup child links
        node->SetupChildLinks(front, back);

        ++ mBvhStats.splits;

        //UpdateViewCellsPvs(front, *frontData.mRays);
        //UpdateViewCellsPvs(back, *backData.mRays);

        ////////////////////////////////////

        //ProcessMultipleRefs(front);
    //ProcessMultipleRefs(back);

        frontData.mNode = front;
        backData.mNode = back;

        back->mObjects = backObjects;
        front->mObjects = frontObjects;

        AssociateObjectsWithLeaf(back);
        AssociateObjectsWithLeaf(front);
   

        // compute probability, i.e., volume of seen view cells
        frontData.mVolume = frontData.mBoundingBox.GetVolume();
        backData.mVolume =  backData.mBoundingBox.GetVolume();

        return node;
}


BvhNode *BvHierarchy::Subdivide(SplitQueue &tQueue,
                                                                SubdivisionCandidate *splitCandidate,
                                                                const bool globalCriteriaMet)
{
        BvhSubdivisionCandidate *sc = dynamic_cast<BvhSubdivisionCandidate *>(splitCandidate);
        BvhTraversalData &tData = sc->mParentData;

        BvhNode *newNode = tData.mNode;

        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
        {       
                BvhTraversalData tFrontData;
                BvhTraversalData tBackData;

                //-- continue subdivision
                
                // create new interior node and two leaf node
                newNode = SubdivideNode(sc->mFrontObjects,
                                                                sc->mBackObjects,
                                                                tData, 
                                                                tFrontData, 
                                                                tBackData);
        
                const int maxCostMisses = sc->mMaxCostMisses;

        // how often was max cost ratio missed in this branch?
                tFrontData.mMaxCostMisses = maxCostMisses;
                tBackData.mMaxCostMisses = maxCostMisses;
                
                mTotalCost -= sc->GetRenderCostDecrease();

                // subdivision statistics
                if (1) EvalSubdivisionStats(*sc);

                //-- push the new split candidates on the queue
                
                BvhSubdivisionCandidate *frontCandidate = new BvhSubdivisionCandidate(tFrontData);
                BvhSubdivisionCandidate *backCandidate = new BvhSubdivisionCandidate(tBackData);

                EvalSubdivisionCandidate(*frontCandidate);
                EvalSubdivisionCandidate(*backCandidate);

                tQueue.Push(frontCandidate);
                tQueue.Push(backCandidate);

                // delete old leaf node
                DEL_PTR(tData.mNode);
        }


        //-- terminate traversal
        if (newNode->IsLeaf())
        {
                EvaluateLeafStats(tData);
        
                const bool mStoreRays= true;

                //-- store additional info
                if (mStoreRays)
                {
                        BvhLeaf *leaf = dynamic_cast<BvhLeaf *>(newNode);
                        CollectRays(leaf->mObjects, leaf->mVssRays);
                }
        }
        
        //-- cleanup
        tData.Clear();
        
        return newNode;
}


void BvHierarchy::EvalSubdivisionCandidate(BvhSubdivisionCandidate &splitCandidate)
{
        ObjectContainer frontObjects, backObjects;

        // compute best object partition
        const float ratio =
                SelectObjectPartition(splitCandidate.mParentData, frontObjects, backObjects);

        const bool success = ratio < mTermMaxCostRatio;

        float oldRenderCost;

        // compute global decrease in render cost
        const float renderCostDecr = EvalRenderCostDecrease(frontObjects, 
                                                                                                                backObjects,
                                                                                                                splitCandidate.mParentData,
                                                                                                                oldRenderCost);

        splitCandidate.SetRenderCostDecrease(renderCostDecr);

#if 0
        const float priority = (float)-data.mDepth;
#else   
        // take render cost of node into account 
        // otherwise danger of being stuck in a local minimum!!
        const float factor = mRenderCostDecreaseWeight;
        const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;

#endif

        // compute global decrease in render cost
        splitCandidate.SetPriority(priority);

        splitCandidate.mMaxCostMisses = 
                success ? splitCandidate.mParentData.mMaxCostMisses : 
        splitCandidate.mParentData.mMaxCostMisses + 1;
}


inline bool BvHierarchy::LocalTerminationCriteriaMet(const BvhTraversalData &data) const
{
        // matt: TODO
        return ( 
                //(data.mNode->mObjects.size() < mTermMinObjects) ||
                //(data.mProbability <= mTermMinProbability) ||
                (data.mDepth >= mTermMaxDepth)
                 );
}


inline bool BvHierarchy::GlobalTerminationCriteriaMet(const BvhTraversalData &data) const
{
        // matt: TODO
        return (
                (mBvhStats.Leaves() >= mTermMaxLeaves)
                //mOutOfMemory || 
                //(mGlobalCostMisses >= mTermGlobalCostMissTolerance)
                );
}


void BvHierarchy::EvaluateLeafStats(const BvhTraversalData &data)
{
        // the node became a leaf -> evaluate stats for leafs
        BvhLeaf *leaf = data.mNode;
        
        ++ mCreatedLeaves;

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

        if (data.mVolume <= mTermMinVolume)
                ++ mBvhStats.minProbabilityNodes;
        
        // accumulate depth to compute average depth
        mBvhStats.accumDepth += data.mDepth;
 
        //      if ((int)(leaf->mObjects.size()) < mTermMinCost)
        //     ++ mOspStats.minCostNodes;
 
        if ((int)(leaf->mObjects.size()) > mBvhStats.maxObjectRefs)
                mBvhStats.maxObjectRefs = (int)leaf->mObjects.size();

#ifdef _DEBUG
        Debug << "BSP stats: "
                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
                  << "Prob: " << data.mProbability << " (min: " << mTermMinProbability << ")\n";
#endif
}


float BvHierarchy::EvalLocalCostHeuristics(const BvhTraversalData &tData,
                                                                                   const int axis,
                                                                                   ObjectContainer &objectsFront,
                                                                                   ObjectContainer &objectsBack)
{
        // the relative cost ratio
        float ratio = 99999999.0f;

        // go through the lists, count the number of objects left and right
        // and evaluate the following cost funcion:
        // C = ct_div_ci  + (ol + or)/queries
        
        const float minBox = tData.mBoundingBox.Min(axis);
        const float maxBox = tData.mBoundingBox.Max(axis);

        const float sizeBox = maxBox - minBox;

        float minBand = minBox + mSplitBorder * (maxBox - minBox);
        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);

        //-- sort so we can use a sweep
        SortSubdivisionCandidates(tData, axis, minBand, maxBand);

        int numViewCells;

        float totalVol = PrepareHeuristics(tData);
        float voll = 0;
        float volr = totalVol;

        /// this is kind of a reverse pvs
        const int totalObjects = (int)tData.mNode->mObjects.size();
        
        int objectsl = 0;
        int objectsr = totalObjects;

        float sum = (float)totalVol * objectsr;


        bool splitPlaneFound = false;

        /////////////////////////////////
        // the parameters for the current optimum

        float volBack = voll;
        float volFront = volr;

        int nObjectsBack = objectsl;
        int nObjectsFront = objectsr;
        
        float minSum = 1e20f;

        /////////////////////////////

        debugVol = 0;
        const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();

        /////////////////////////////
        // the sweep heuristics

        // mail the objects on the left side
        Intersectable::NewMail();
        
        //-- traverse through events and find best split plane
        vector<SortableEntry>::const_iterator ci, ci_end = mSubdivisionCandidates->end();

        for (ci = mSubdivisionCandidates->begin(); ci != ci_end; ++ ci)
        {
                Intersectable *object = (*ci).mObject;

                EvalHeuristicsContribution(tData.mNode, *ci, voll, volr, objectsl);
                objectsr = totalObjects - objectsl;

                // Note: sufficient to compare size of bounding boxes of front and back side?

                if (((*ci).mPos >= minBand) && ((*ci).mPos <= maxBand))
                {
                        // voll = view cells that see only left node (i.e., left pvs)
                        // volr = view cells that see only right node (i.e., right pvs)
                        // rest = view cells that see both nodes (i.e., total pvs)
            sum = voll * objectsl + volr * objectsr + (totalVol - voll - volr) * totalObjects;

                        float currentPos;
                        
                        // HACK: current positition is BETWEEN visibility events
                        if (1 && ((ci + 1) != ci_end))
                                currentPos = ((*ci).mPos + (*(ci + 1)).mPos) * 0.5f;
                        else
                                currentPos = (*ci).mPos;                        

                        if ((totalVol - voll - volr - debugVol) / viewSpaceVol > 0.0001)
                                Debug << "front and back volume: " << (totalVol - voll - volr) / viewSpaceVol << " error: " << (totalVol - voll - volr - debugVol) / viewSpaceVol << endl;
#if 0                           
                        Debug << "pos: " << currentPos 
                                 << "\t (pvsl: " << pvsl << ", pvsr: " << pvsr << ")"
                                 << "\t (voll: " << voll << ", volr: " << volr << ")"
                                 << "\t (vcl: " << vcl << ", vcr: " << vcr << ", nvc: " << numViewCells << ")"
                                 << "\t sum: " << sum << endl;
                                
#endif
                        if (sum < minSum)
                        {
                                splitPlaneFound = true;

                                minSum = sum;

                                nObjectsBack = objectsl;
                                nObjectsFront = objectsr;
                                
                                volBack = voll;
                                volFront = volr;
                        }
                }
        }
        
        //-- compute cost
        const float frontAndBackVol = totalVol - volFront - volBack;

        const float oldRenderCost = (float)totalObjects * totalVol + Limits::Small;
        const float newRenderCost = (float)nObjectsFront * volFront + 
                                                                (float)nObjectsBack * volBack +
                                                                (float)totalObjects * frontAndBackVol;

        if (splitPlaneFound)
        {
                ratio = newRenderCost / oldRenderCost;
        }

        //Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" 
        //        << position << " t=" << (position - minBox) / (maxBox - minBox)
        //      << "\t pb=(" << volBack << ")\t pf=(" << volFront << ")" << endl;

        Debug << "\n§§§§ eval local cost §§§§" << endl
                  << "back pvs: " << nObjectsBack << " front pvs: " << nObjectsFront << " total pvs: " << totalObjects << endl 
                  << "back p: " << volBack / viewSpaceVol << " front p " << volFront / viewSpaceVol << " p: " << totalVol / viewSpaceVol << endl
                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
                  << "render cost decrease: " << oldRenderCost / viewSpaceVol - newRenderCost / viewSpaceVol << endl;

        if (oldRenderCost < newRenderCost)
                Debug << "\nwarning!!:\n" << "old rc: " << oldRenderCost * viewSpaceVol << " new rc: " << newRenderCost * viewSpaceVol << endl;

        // assign objects
        ObjectContainer::const_iterator oit, oit_end = tData.mNode->mObjects.end();

        for (oit = tData.mNode->mObjects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;

                if (obj->Mailed()) // belongs to back objects
                        objectsBack.push_back(obj);
                else
                        objectsFront.push_back(obj);
        }

        return ratio;
}


void BvHierarchy::SortSubdivisionCandidates(const BvhTraversalData &tData,
                                                                          const int axis, 
                                                                          float minBand, 
                                                                          float maxBand)

{
        mSubdivisionCandidates->clear();

        //RayInfoContainer *rays = tData.mRays;
        BvhLeaf *leaf = tData.mNode;

        // compute requested size
        int requestedSize = 0;
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
        for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
                requestedSize += 2 + (int)(*oit)->mVssRays.size() * 2;
        
        // creates a sorted split candidates array
        if (mSubdivisionCandidates->capacity() > 500000 &&
                requestedSize < (int)(mSubdivisionCandidates->capacity() / 10) )
        {
        delete mSubdivisionCandidates;
                mSubdivisionCandidates = new vector<SortableEntry>;
        }

        mSubdivisionCandidates->reserve(requestedSize);

        float pos;

        //-- insert object queries
        //-- These queries can induce a change in pvs size.
        //-- Note that they cannot induce a change in view cell volume, as
        //-- there is no ray associated with these events!!

        for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
        {
                Intersectable *obj = *oit;
                
                Intersectable *object = *oit;
                const AxisAlignedBox3 box = object->GetBox();

                mSubdivisionCandidates->push_back(
                        SortableEntry(SortableEntry::OBJECT,
                                                  box.Min(axis),
                                                  object,
                                                  NULL)
                                                  );


                //-- insert ray queries
                //-- we are intersested only in rays which intersect an object that
                //-- is part of the kd node because they can induce a change in view cell
                //-- volume on the left and the right part. 
                //-- Note that they cannot induce a change in pvs size!!
                
                VssRayContainer::const_iterator rit, rit_end = obj->mVssRays.end();

                for (rit = obj->mVssRays.begin(); rit < rit_end; ++ rit)
                {
                        VssRay *ray = (*rit);
                        const bool positive = ray->HasPosDir(axis);
        
                        pos = ray->mTermination[axis];

                        mSubdivisionCandidates->push_back(
                                SortableEntry(SortableEntry::VIEWCELLS, 
                                pos, 
                                ray->mTerminationObject, 
                                ray)
                                );
                }
        }

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


const BvhStatistics &BvHierarchy::GetStatistics() const
{
        return mBvhStats;
}


void BvHierarchy::EvalRayContribution(const VssRay &ray, 
                                                                          float &volLeft, 
                                                                          float &volRight)
{
        ViewCellContainer viewCells;
        mVspTree->GetViewCells(ray, viewCells);

        /// classify view cells and compute volume contri accordingly
        /// possible view cell classifications:
        /// view cell mailed => view cell can be seen from left child node
        /// view cell counter > 0 view cell can be seen from right child node
        /// combined: view cell volume belongs to both nodes
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
        
        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
        {
                // view cells can also be seen from left child node
                ViewCell *viewCell = *vit;

                const float vol = viewCell->GetVolume();

                if (!viewCell->Mailed())
                {
                        viewCell->Mail();

                        // we now see view cell from both nodes 
                        // => remove ref to right child node
                        volRight -= vol;
                        debugVol += vol;
                }

                // last reference into the right node
                if (-- viewCell->mCounter == 0)
                {
                        // view cell was previously seen from both nodes  =>
                        // contributes only to left node now
                        volLeft += vol;
                        debugVol -= vol;
                }
        }
}


float BvHierarchy::PrepareHeuristics(const VssRay &ray)
{
        float vol = 0;

        ViewCellContainer viewCells;
        mVspTree->GetViewCells(ray, viewCells);
        
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = (*vit);

                if (!vc->Mailed())
                {
                        //Debug << "single vol: "<< vc->GetVolume() << endl;
                        vc->Mail();
                        vc->mCounter = 0;
                        vol += vc->GetVolume();
                }
                
                // view cell volume already added => just increase counter
                ++ vc->mCounter;
        }

        return vol;
}


float BvHierarchy::PrepareHeuristics(const BvhTraversalData &tData)
{       
        float vol = 0;
        debugVol = 0;
        
        ViewCell::NewMail();

        BvhLeaf *leaf = tData.mNode;

        VssRayContainer rays;
        CollectRays(leaf->mObjects, rays);

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

        for (rit = rays.begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit);

                const float volContri = PrepareHeuristics(*ray);

                // if hitpoint with one of the objects is inside this node, we
                // evaluate the volume of the view cells seen by this ray
                if (ray->mTerminationObject)
                {
            vol += volContri;
                }

                // count double if both hit points are within the kd node
                if (0 && ray->mOriginObject)
                {
                        vol += volContri;
                }
        }

        return vol;
}
///////////////////////////////////////////////////////////


void BvHierarchy::EvalHeuristicsContribution(BvhLeaf *leaf,
                                                                                         const SortableEntry &ci,
                                                                                         float &volLeft, 
                                                                                         float &volRight, 
                                                                                         int &objectsLeft)
{
        Intersectable *obj = ci.mObject;
        VssRay *ray = ci.mRay;

        // switch between different types of events
        switch (ci.mType) 
        {
                case SortableEntry::OBJECT:
                        cout << "o";
                        ci.mObject->Mail();
                        ++ objectsLeft;
                        break;

                // compute volume contribution from view cells
                case SortableEntry::VIEWCELLS:
                        cout << "v";
                        // process ray if the hit point with termination / origin object 
                        // is inside this kd leaf
                        EvalRayContribution(*ray, volLeft, volRight);
                        break;

                default:
                        cout << "should not come here" << endl;
                        break;
        }
        //cout << "vol left: " << volLeft << " vol right " << volRight << endl;
}


void BvHierarchy::SetViewCellsManager(ViewCellsManager *vcm)
{
        mViewCellsManager = vcm;
}


AxisAlignedBox3 BvHierarchy::GetBoundingBox() const
{
        return mBoundingBox;
}


float BvHierarchy::SelectObjectPartition(const BvhTraversalData &tData,
                                                                                 ObjectContainer &frontObjects,
                                                                                 ObjectContainer &backObjects)
{
        ObjectContainer nFrontObjects[3];
        ObjectContainer nBackObjects[3];

        float nCostRatio[3];

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

        if (mOnlyDrivingAxis)
        {
                sAxis = box.Size().DrivingAxis();
        }
        
        // -- evaluate split cost for all three axis
        
        for (int axis = 0; axis < 3; ++ axis)
        {
                if (!mOnlyDrivingAxis || (axis == sAxis))
                {
                        if (1 || mUseCostHeuristics)
                        {
                                //-- partition objects using heuristics
                                
                                nCostRatio[axis] =
                                        EvalLocalCostHeuristics(
                                                                                        tData,
                                                                                        axis,
                                                                                        nFrontObjects[axis],
                                                                                        nBackObjects[axis]);                    
                        }
                        if (bestAxis == -1)
                        {
                                bestAxis = axis;
                        }
                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
                        {
                                bestAxis = axis;
                        } 
                }
        }

        //-- assign values
        
        frontObjects = nFrontObjects[bestAxis];
        backObjects = nFrontObjects[bestAxis];

        Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
        return nCostRatio[bestAxis];
}


void BvHierarchy::AssociateObjectsWithRays(const VssRayContainer &rays)
{

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

    for (rit = rays.begin(); rit != rays.end(); ++ rit)
        {
                VssRay *ray = (*rit);

                if (ray->mTerminationObject)
                {
                        ray->mTerminationObject->mVssRays.push_back(ray);
                }

                if (0 && ray->mOriginObject)
                {
                        ray->mOriginObject->mVssRays.push_back(ray);
                }
        }
}


void BvHierarchy::EvalSubdivisionStats(const SubdivisionCandidate &sc)
{
        const float costDecr = sc.GetRenderCostDecrease();

        AddSubdivisionStats(mBvhStats.Leaves(),
                                                costDecr,
                                                mTotalCost
                                                );
}


void BvHierarchy::AddSubdivisionStats(const int leaves,
                                                                          const float renderCostDecr,
                                                                          const float totalRenderCost)
{
        mSubdivisionStats 
                        << "#Leaves\n" << leaves << endl
                        << "#RenderCostDecrease\n" << renderCostDecr << endl 
                        << "#TotalRenderCost\n" << totalRenderCost << endl;
                        //<< "#AvgRenderCost\n" << avgRenderCost << endl;
}


void BvHierarchy::CollectRays(const ObjectContainer &objects,
                                                          VssRayContainer &rays) const
{
        VssRay::NewMail();
        
        ObjectContainer::const_iterator oit, oit_end = objects.end();

        // evaluate reverse pvs and view cell volume on left and right cell
        // note: should I take all leaf objects or rather the objects hit by rays?
        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;
                
                VssRayContainer::const_iterator rit, rit_end = obj->mVssRays.end();

                for (rit = obj->mVssRays.begin(); rit < rit_end; ++ rit)
                {
                        VssRay *ray = (*rit);

                        if (!ray->Mailed())
                        {
                                ray->Mail();
                                rays.push_back(ray);
                        }
                }
        }
}


void BvHierarchy::MailViewCells(const ObjectContainer &objects, 
                                                                const bool isFront,
                                                                ViewCellContainer &collectedViewCells) const
{
        VssRayContainer rays;
        CollectRays(objects, rays);
        
        VssRayContainer::const_iterator rit, rit_end = rays.end();

        for (rit = rays.begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit);

                // view cell is assigned to left and / or right child
                ViewCellContainer viewCells;
                mVspTree->GetViewCells(*ray, viewCells);
                        
                ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

                // mail view cell
                for (vit = viewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;

                        if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
                                collectedViewCells.push_back(vc);

                        if (isFront)
                        {
                                if (!vc->Mailed())
                                        vc->Mail();
                        }
                        else
                        {
                                if (!vc->Mailed())
                                        vc->Mail(1);
                                else
                                        vc->Mail(2);
                        }
                }                               
        }
}


float BvHierarchy::EvalRenderCostDecrease(const ObjectContainer &objectsFront,
                                                                                  const ObjectContainer &objectsBack,
                                                                                  const BvhTraversalData &tData,
                                                                                  float &normalizedOldRenderCost) const
{
        // probability that view point lies in back / front node
        float pOverall = 0;
        float pFront = 0;
        float pBack = 0;
        float pFrontAndBack = 0;

        const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();

        BvhLeaf *leaf = tData.mNode;
        
        // sum up volume seen from the objects of left and right children
        // => the volume is the weight for the render cost equation
        ViewCell::NewMail(3);

        ViewCellContainer collectedViewCells;
        MailViewCells(objectsFront, true, collectedViewCells);
        MailViewCells(objectsBack, false, collectedViewCells);

        ViewCellContainer::const_iterator vit, vit_end = collectedViewCells.end();

        // evaluate view cells volume contribution with respect to the mail id
        for (vit = collectedViewCells.begin(); vit != vit_end; ++ vit)
        {
                AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
        }

        const int totalObjects = (int)leaf->mObjects.size();
        const int nObjectsFront = (int)objectsFront.size();
        const int nObjectsBack = (int)objectsBack.size();

        // these are non-overlapping sets
        pOverall = pFront + pBack + pFrontAndBack;

        //-- pvs rendering heuristics
        const float oldRenderCost = pOverall * totalObjects;
        const float newRenderCost = nObjectsFront * pFront + nObjectsBack * pBack + totalObjects * pFrontAndBack;

        // normalize volume with view space volume
        const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;

        Debug << "\n(((( eval render cost decrease ))))" << endl
                  << "back objects: " << nObjectsBack << " front objects " << nObjectsFront << " total objects: " << totalObjects << endl 
                  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol 
                  << " front and back p " << pFrontAndBack / viewSpaceVol << " p: " << pOverall / viewSpaceVol << endl
                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
                  << "render cost decrease: " << renderCostDecrease << endl;

        normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

        if (oldRenderCost < newRenderCost * 0.99)
                Debug << "\nwarning2!!:\n" << "old rc: " << oldRenderCost * viewSpaceVol << " new rc: " << newRenderCost * viewSpaceVol << endl;
        
        return renderCostDecrease;
}


AxisAlignedBox3 BvHierarchy::ComputeBoundingBox(const ObjectContainer &objects)
{
        AxisAlignedBox3 box;
        box.Initialize();

        ObjectContainer::const_iterator oit, oit_end = objects.end();

        //-- compute bounding box
        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;

                // compute bounding box of view space
                mBoundingBox.Include(obj->GetBox());
        }
        return box;
        
}


void BvHierarchy::CollectLeaves(vector<BvhLeaf *> &leaves) const
{
        stack<BvhNode *> nodeStack;
        nodeStack.push(mRoot);

        while (!nodeStack.empty()) 
        {
                BvhNode *node = nodeStack.top();
                nodeStack.pop();
                if (node->IsLeaf()) 
                {
                        BvhLeaf *leaf = (BvhLeaf *)node;
                        leaves.push_back(leaf);
                }
                else 
                {
                        BvhInterior *interior = (BvhInterior *)node;

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


AxisAlignedBox3 BvHierarchy::GetBoundingBox(BvhNode *node) const
{
        return node->GetBoundingBox();
}


void BvHierarchy::CollectViewCells(BvhLeaf *leaf, ViewCellContainer &viewCells)
{
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

        ViewCell::NewMail();

        // loop through all object and collect view cell pvs of this node
        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;

                ViewCellPvsMap::const_iterator vit, vit_end = obj->mViewCellPvs.mEntries.end();

                for (vit = obj->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
                {
            ViewCell *vc = (*vit).first;

                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                viewCells.push_back(vc);
                        }
                }
        }
}


void BvHierarchy::CollectDirtyCandidates(BvhSubdivisionCandidate *sc, 
                                                                                 vector<SubdivisionCandidate *> &dirtyList)
{
        BvhTraversalData &tData = sc->mParentData;
        BvhLeaf *node = tData.mNode;
        
        ViewCell::NewMail();

        ViewCellContainer viewCells;
        ViewCellContainer tmpViewCells;

        VssRayContainer rays;
        CollectRays(node->mObjects, rays);

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

        // find all view cells associated with the rays, add them to view cells
        for (rit = rays.begin(); rit != rit_end; ++ rit)
        {
                VssRay *ray = (*rit);
                mVspTree->GetViewCells(*ray, tmpViewCells);

                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();

                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
                {
                        VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);

                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                viewCells.push_back(vc);
                        }
                }
        }


        // split candidates holding this view cells 
        // are thrown into dirty list
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
        {
                VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);

                VspLeaf *leaf = vc->mLeaf;
                dirtyList.push_back(leaf->GetSubdivisionCandidate());
        }
}


BvhNode *BvHierarchy::GetRoot() const
{
        return mRoot;
}


bool BvHierarchy::IsObjectInLeaf(BvhLeaf *leaf, Intersectable *object) const
{
        ObjectContainer::const_iterator oit =
                lower_bound(leaf->mObjects.begin(), leaf->mObjects.end(), object, ilt);
                                
        // objects sorted by id
        if ((oit != leaf->mObjects.end()) && ((*oit)->GetId() == object->GetId()))
        {
                return true;
        }
        else
        {
                return false;
        }
}


BvhLeaf *BvHierarchy::GetLeaf(Intersectable *object, BvhNode *node) const
{
        // rather use the simple version
        return object->mBvhLeaf;

        ///////////////////////////////////////
        // start from root of tree
        if (node == NULL)
        {
                node = mRoot;
        }

        vector<BvhLeaf *> leaves;

        stack<BvhNode *> nodeStack;
        nodeStack.push(node);
  
        BvhLeaf *leaf = NULL;
  
        while (!nodeStack.empty())  
        {
                BvhNode *node = nodeStack.top();
                nodeStack.pop();
        
                if (node->IsLeaf()) 
                {
                        leaf = dynamic_cast<BvhLeaf *>(node);

                        if (IsObjectInLeaf(leaf, object))
                                return leaf;
                } 
                else    
                {       
                        // find point
                        BvhInterior *interior = dynamic_cast<BvhInterior *>(node);
        
                        if (interior->GetBack()->GetBoundingBox().Includes(object->GetBox()))
                        {
                                nodeStack.push(interior->GetBack());
                        }
                        
                        // search both sides as we are using bounding volumes
                        if (interior->GetFront()->GetBoundingBox().Includes(object->GetBox()))
                        {
                                nodeStack.push(interior->GetFront());
                        }
                }
        }
  
        return leaf;
}


BvhIntersectable *BvHierarchy::GetOrCreateBvhIntersectable(BvhNode *node)
{
        // search nodes
        std::map<BvhNode *, BvhIntersectable *>::
                const_iterator it = mBvhIntersectables.find(node);

        if (it != mBvhIntersectables.end()) 
        {
                return (*it).second;
        }

        // not in map => create new entry
        BvhIntersectable *bvhObj = new BvhIntersectable(node);
        mBvhIntersectables[node] = bvhObj;

        return bvhObj;
}


void BvHierarchy::AddViewCellVolumeContri(ViewCell *vc,
                                                                                  float &frontVol,
                                                                                  float &backVol,
                                                                                  float &frontAndBackVol) const
{
        if (vc->Mailed())
        {
                frontVol += vc->GetVolume();
        }
        else if (vc->Mailed(1))
        {
                backVol += vc->GetVolume();
        }
        else if (vc->Mailed(2))
        {
                frontAndBackVol += vc->GetVolume();
        }
}

/*
int BvHierarchy::UpdateViewCellsPvs(BvhLeaf *leaf, 
                                                                const RayInfoContainer &rays) const

{
        MailablePvsData::NewMail();
        
        ViewCellContainer touchedViewCells;
        CollectTouchedViewCells(rays, touchedViewCells);

        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

        for (oit = leaf->mObjects.begin(); oit < oit_end; ++ oit)
        {
                Intersectable *obj = *oit;
                ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

                // traverse through view cells and classify them according
                // to them being seen from to back / front / front and back node
                for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;
                        float contri;
                        AddViewCellToObjectPvs(obj, vc, contri, true);
                }
        }

        return 0;
}


int BvHierarchy::RemoveParentViewCellsPvs(BvhLeaf *leaf, 
                                                                          const RayInfoContainer &rays
                                                                          ) const

{
        MailablePvsData::NewMail();
        
        ViewCellContainer touchedViewCells;
        CollectTouchedViewCells(rays, touchedViewCells);

        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;

                // traverse through view cells and classify them according
                // to them being seen from to back / front / front and back node
        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

                for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;
                        
                        MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

                        if (vdata && !vdata->Mailed())
                        {
                                vdata->Mail();
                                obj->mViewCellPvs.RemoveSample(vc, 1);
                        }
                }
        }

        return 0;
}
*/

bool BvHierarchy::Export(OUT_STREAM &stream)
{
        ExportNode(mRoot, stream);

        return true;
}


void BvHierarchy::ExportNode(BvhNode *node, OUT_STREAM &stream)
{
        if (node->IsLeaf())
        {
                BvhLeaf *leaf = dynamic_cast<BvhLeaf *>(node);

                stream << "<Leaf ";
                stream << "objects=\"";
                
                //-- export objects in kd leaves
                //if (mExportObjects) ExportObjects(leaf, stream);
                
                stream << "\" />" << endl;
                //stream << " </Leaf>" << endl;
        }
        else
        {       
                BvhInterior *interior = dynamic_cast<BvhInterior *>(node);
        
                stream << "<Interior box=\"" << interior->GetBoundingBox() << "\">" << endl;

                ExportNode(interior->GetBack(), stream);
                ExportNode(interior->GetFront(), stream);

                stream << "</Interior>" << endl;
        }
}


float BvHierarchy::EvalViewCellsVolume(BvhLeaf *leaf) const
{
        float vol = 0;

        VssRayContainer rays;
        CollectRays(leaf->mObjects, rays);

        ViewCell::NewMail();

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

        for (rit = rays.begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit);
        
                ViewCellContainer viewCells;
                mVspTree->GetViewCells(*ray, viewCells);
        
                ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

                for (vit = viewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = (*vit);

                        if (!vc->Mailed())
                        {
                                //Debug << "single vol: "<< vc->GetVolume() << endl;
                                vc->Mail();
                                vol += vc->GetVolume();
                        }
                }
        }

        return vol;
}


SubdivisionCandidate *BvHierarchy::PrepareConstruction(const VssRayContainer &sampleRays,
                                                                                                           ObjectContainer &objects,
                                                                                                           AxisAlignedBox3 *forcedObjectSpace
                                                                                                           //, RayInfoContainer &rays
                                                                                                           )
{
        // store pointer to this tree
        BvhSubdivisionCandidate::sBvHierarchy = this;
        mBvhStats.nodes = 1;

        // compute bounding box from objects
        if (forcedObjectSpace)
                mBoundingBox = *forcedObjectSpace;
        else
                mBoundingBox = ComputeBoundingBox(objects);

        mTermMinVolume *= mBoundingBox.GetVolume();
        mGlobalCostMisses = 0;

        // sort objects by id in order to have sorted objects in
        // the leaf nodes
        stable_sort(objects.begin(), objects.end(), ilt);

        //-- create new root

        BvhLeaf *bvhleaf = new BvhLeaf(mBoundingBox, NULL, (int)objects.size());
        bvhleaf->mObjects = objects;
        mRoot = bvhleaf;

        // only rays intersecting objects in node are interesting
        AssociateObjectsWithRays(sampleRays);
        // associate root with current objects
        AssociateObjectsWithLeaf(bvhleaf);

        //-- add first candidate for object space partition

        // create osp traversal data
        BvhTraversalData oData(bvhleaf, 0, mBoundingBox.GetVolume());

        //-- first split candidate
        BvhSubdivisionCandidate *oSubdivisionCandidate = 
                new BvhSubdivisionCandidate(oData);

        //UpdateViewCellsPvs(kdleaf, rays);

        EvalSubdivisionCandidate(*oSubdivisionCandidate);

// probabilty is voume of all "seen" view cells
#if 1
        const float prop = EvalViewCellsVolume(bvhleaf);
#else
        const float prop = GetBoundingBox().GetVolume();
#endif
        mTotalCost = (float)objects.size() * prop / 
                mVspTree->GetBoundingBox().GetVolume();

        EvalSubdivisionStats(*oSubdivisionCandidate);

        return oSubdivisionCandidate;
}



}