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

#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"
#include "HierarchyManager.h"


namespace GtpVisibilityPreprocessor {


#define PROBABILIY_IS_BV_VOLUME 1
#define USE_FIXEDPOINT_T 0
#define USE_VOLUMES_FOR_HEURISTICS 1
#define TEST_POWERPLANT 0


BvHierarchy *BvHierarchy::BvhSubdivisionCandidate::sBvHierarchy = NULL;


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


/// sorting operator
inline static bool smallerSize(Intersectable *obj1, Intersectable *obj2)
{
        return obj1->GetBox().SurfaceArea() < obj2->GetBox().SurfaceArea();
}



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


BvhNode::BvhNode(): 
mParent(NULL), 
mTimeStamp(0),
mRenderCost(-1)
{
}

BvhNode::BvhNode(const AxisAlignedBox3 &bbox):
mParent(NULL), 
mBoundingBox(bbox), 
mTimeStamp(0),
mRenderCost(-1)
{
}


BvhNode::BvhNode(const AxisAlignedBox3 &bbox, BvhInterior *parent):
mBoundingBox(bbox), 
mParent(parent), 
mTimeStamp(0),
mRenderCost(-1)
{
}


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


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


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


int BvhNode::GetRandomEdgePoint(Vector3 &point,
                                                                Vector3 &normal)
{
        // get random edge
        const int idx = Random(12); 
        Vector3 a, b;
        mBoundingBox.GetEdge(idx, &a, &b);
        
        const float w = RandomValue(0.0f, 1.0f);

        point = a * w + b * (1.0f - w);

        // TODO
        normal = Vector3(0);

        return idx;
}



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


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox):
BvhNode(bbox), 
mSubdivisionCandidate(NULL),
mGlList(0)
{
  mActiveNode = this;
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent):
  BvhNode(bbox, parent),
  mGlList(0)
  
{
        mActiveNode = this;
}


BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, 
                                 BvhInterior *parent, 
                                 const int numObjects):
  BvhNode(bbox, parent),
  mGlList(0)

{
        mObjects.reserve(numObjects);
        mActiveNode = this;
}


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


BvhLeaf::~BvhLeaf() 
{ 
}


void BvhLeaf::CollectObjects(ObjectContainer &objects)
{
        ObjectContainer::const_iterator oit, oit_end = mObjects.end();
        for (oit = mObjects.begin(); oit != oit_end; ++ oit)
        {
                objects.push_back(*oit);
        }
}



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


BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox):
BvhNode(bbox), mFront(NULL), mBack(NULL)
{
}


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


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;
}


void BvhInterior::CollectObjects(ObjectContainer &objects)
{
        mFront->CollectObjects(objects);
        mBack->CollectObjects(objects);
}


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


BvHierarchy::BvHierarchy():
mRoot(NULL),
mTimeStamp(1),
mIsInitialSubdivision(false)
{
        ReadEnvironment();
        mSubdivisionCandidates = new SortableEntryContainer;
//      for (int i = 0; i < 4; ++ i)
//              mSortedObjects[i] = NULL;
}


BvHierarchy::~BvHierarchy()
{
        // delete the local subdivision candidates
        DEL_PTR(mSubdivisionCandidates);

        // delete the presorted objects
        /*for (int i = 0; i < 4; ++ i)
        {
                DEL_PTR(mSortedObjects[i]);
        }*/
        
        // delete the tree
        DEL_PTR(mRoot);
}


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

        //////////////////////////////
        //-- termination criteria for autopartition

        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.maxDepth", mTermMaxDepth);
        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.maxLeaves", mTermMaxLeaves);
        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.minObjects", mTermMinObjects);
        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.minRays", mTermMinRays);
        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.minProbability", mTermMinProbability);
    Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.missTolerance", mTermMissTolerance);


        //////////////////////////////
        //-- max cost ratio for early tree termination

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


        //////////////////////////////
        //-- factors for subdivision heuristics

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

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

        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
        Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.useGlobalSorting", mUseGlobalSorting);
        Environment::GetSingleton()->GetIntValue("BvHierarchy.minRaysForVisibility", mMinRaysForVisibility);
        Environment::GetSingleton()->GetIntValue("BvHierarchy.maxTests", mMaxTests);
        Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.useInitialSubdivision", mApplyInitialPartition);
        Environment::GetSingleton()->GetIntValue("BvHierarchy.Construction.Initial.minObjects", mInitialMinObjects);
        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.Initial.maxAreaRatio", mInitialMaxAreaRatio);
        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.Initial.minArea", mInitialMinArea);

        //mMemoryConst = (float)(sizeof(VspLeaf) + sizeof(VspViewCell));
        //mMemoryConst = (float)sizeof(BvhLeaf);
        mMemoryConst = 16;//(float)sizeof(ObjectContainer);

    mUseBboxAreaForSah = true;

        /////////////
        //-- debug output

        Debug << "******* Bvh hierarchy options ******** " << endl;
    Debug << "max depth: " << mTermMaxDepth << endl;
        Debug << "min probabiliy: " << mTermMinProbability<< 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 << "use surface area heuristics: " << mUseSah << endl;
        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
        //Debug << "split borders: " << mSplitBorder << endl;
        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
        Debug << "use global sort: " << mUseGlobalSorting << endl;
        Debug << "minimal rays for visibility: " << mMinRaysForVisibility << endl;
        Debug << "bvh mem const: " << mMemoryConst << endl;
        Debug << "apply initial partition: " << mApplyInitialPartition << endl;
        Debug << "min objects: " << mInitialMinObjects << endl;
        Debug << "max area ratio: " << mInitialMaxAreaRatio << endl;
        Debug << "min area: " << mInitialMinArea << endl;

        Debug << endl;
}


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

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


static int CountRays(const ObjectContainer &objects)
{
        int nRays = 0;

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

        // warning: not exact number (there can be rays counted twice)
        // otherwise we would have to traverse trough all rays
        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                nRays += (int)(*oit)->GetOrCreateRays()->size();
        }

        return nRays;
}

                                                                        
float BvHierarchy::GetViewSpaceVolume() const
{
        return mViewCellsManager->GetViewSpaceBox().GetVolume();
}


void BvHierarchy::UpdateViewCells(const BvhSubdivisionCandidate &sc)
{
        ViewCellContainer viewCells, frontViewCells, backViewCells;
        
        CollectViewCells(*sc.mParentData.mSampledObjects, viewCells, false, false);
        CollectViewCells(sc.mSampledFrontObjects, frontViewCells, false, false);
        CollectViewCells(sc.mSampledBackObjects, backViewCells, false, false);

        const int frontTri = (int)sc.mFrontObjects.size();
        const int backTri = (int)sc.mBackObjects.size();
        const int totalTri = (int)(*sc.mParentData.mSortedObjects[0]).size();

        //cout << "totalTri: " << totalTri << " f: " << frontTri << " back: " << backTri << endl;

        ViewCell::NewMail(3);

        // mail view cells which can see front object
        ViewCellContainer::const_iterator fit, fit_end = frontViewCells.end();

        for (fit = frontViewCells.begin(); fit != fit_end; ++ fit)
        {
                (*fit)->Mail(0);
        }

        // mail view cells which can see back or both objects
        ViewCellContainer::const_iterator bit, bit_end = backViewCells.end();

        for (bit = backViewCells.begin(); bit != bit_end; ++ bit)
        {
                ViewCell *vc = *bit;

                if (vc->Mailed(0))
                        vc->Mail(2);
                else
                        vc->Mail(1);
        }

        // traverse through view cells and compute changes
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
        
        for (vit = viewCells.begin(); vit != viewCells.end(); ++ vit)
        {
                ViewCell *vc = *vit;

                int vcTri;
                int vcObj;

                int oldVcTri = (int)vc->GetTrianglesInPvs();
                int oldVcObj = vc->GetEntriesInPvs();

                // both objects seen from view cell
                // => no reduction, but an additional pvs entry
                if (vc->Mailed(2))
                {
                        vcTri = oldVcTri;
                        vcObj = oldVcObj + 1;   
                }
                // only back object seen from view cell 
                // => reduction in triangles
                else if (vc->Mailed(1))
                {
                        vcTri = oldVcTri + backTri - totalTri; 
                        vcObj = oldVcObj;   
                }
                else // front object
                {
                        vcTri = oldVcTri + frontTri - totalTri; 
                        vcObj = oldVcObj;
                }

                vc->SetTrianglesInPvs((float)vcTri);
                vc->SetEntriesInPvs(vcObj);

                //cout << "old pvs tri: " << oldVcTri << " new: " << vcTri << endl;
                //cout << "old pvs obj: " << oldVcObj << " new: " << vcObj << endl;
        }
}


void BvHierarchy::TestEvaluation(const BvhSubdivisionCandidate &sc)
{
        ViewCellContainer viewCells, frontViewCells, backViewCells;
        
        CollectViewCells(*sc.mParentData.mSampledObjects, viewCells, false, false);

        // traverse through view cells and compute changes
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
        
        for (vit = viewCells.begin(); vit != viewCells.end(); ++ vit)
        {
                ViewCell *vc = *vit;

                int oldVcTri = (int)vc->GetTrianglesInPvs();
                int oldVcObj = vc->GetEntriesInPvs();

                int nTriangles = 0;
                int nObjects = 0;

                Intersectable::NewMail();

                VspViewCell *vspVc = static_cast<VspViewCell *>(vc);
                VssRayContainer::const_iterator vit, vit_end = vspVc->mLeaves[0]->mVssRays.end();

                for (vit = vspVc->mLeaves[0]->mVssRays.begin(); vit != vit_end; ++ vit)
                {
                        VssRay *ray = *vit;
                        
                        BvhLeaf *obj = ray->mTerminationObject->mBvhLeaf;
                
                        if (!obj->Mailed())
                        {
                                obj->Mail();

                                nTriangles += (int)obj->mObjects.size();
                                nObjects ++;
                        }

                        if (ray->mOriginObject)
                        {
                                obj = ray->mOriginObject->mBvhLeaf;
                        
                                if (!obj->Mailed())
                                {
                                        obj->Mail();
                                        nTriangles += (int)obj->mObjects.size();
                                        nObjects ++;
                                }
                        }
                }

                cout << "old pvs tri: " << oldVcTri << " real: " << nTriangles << endl;
                cout << "old pvs obj: " << oldVcObj << " real: " << nObjects << endl;
        }
}



BvhInterior *BvHierarchy::SubdivideNode(const BvhSubdivisionCandidate &sc,
                                                                                BvhTraversalData &frontData,
                                                                                BvhTraversalData &backData)
{ 
        //TestEvaluation(sc);

#ifdef USE_PERFTIMER   
        // fill view cells cache
        mNodeTimer.Entry();
#endif  

        const BvhTraversalData &tData = sc.mParentData;
        BvhLeaf *leaf = tData.mNode;

        AxisAlignedBox3 parentBox = leaf->GetBoundingBox();

        // update stats: we have two new leaves
        mBvhStats.nodes += 2;

        if (tData.mDepth > mBvhStats.maxDepth)
        {
                mBvhStats.maxDepth = tData.mDepth;
        }

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


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

        AxisAlignedBox3 fbox = EvalBoundingBox(sc.mFrontObjects, &parentBox);
        AxisAlignedBox3 bbox = EvalBoundingBox(sc.mBackObjects, &parentBox);

        BvhLeaf *back = new BvhLeaf(bbox, node, (int)sc.mBackObjects.size());
        BvhLeaf *front = new BvhLeaf(fbox, node, (int)sc.mFrontObjects.size());

        BvhInterior *parent = leaf->GetParent();

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

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

        ++ mBvhStats.splits;

 
        ////////////////////////////////
        //-- fill front and back traversal data with the new values

        frontData.mDepth = backData.mDepth = tData.mDepth + 1;

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

        backData.mSampledObjects = new ObjectContainer();
        frontData.mSampledObjects = new ObjectContainer();

        *backData.mSampledObjects = sc.mSampledBackObjects;
        *frontData.mSampledObjects = sc.mSampledFrontObjects;

        back->mObjects = sc.mBackObjects;
        front->mObjects = sc.mFrontObjects;

        // if the number of rays is too low, no assumptions can be made
        // (=> switch to surface area heuristics?)
        frontData.mNumRays = CountRays(sc.mSampledFrontObjects);
        backData.mNumRays = CountRays(sc.mSampledBackObjects);

        AssociateObjectsWithLeaf(back);
        AssociateObjectsWithLeaf(front);
  
        ////////////
        //-- compute pvs correction to cope with undersampling

        frontData.mPvs = (float)sc.mNumFrontViewCells;
        backData.mPvs = (float)sc.mNumBackViewCells;

        frontData.mCorrectedPvs = sc.mCorrectedFrontPvs;
        backData.mCorrectedPvs = sc.mCorrectedBackPvs;


        // compute probability of this node being visible, 
        // i.e., volume of the view cells that can see this node
        frontData.mVolume = sc.mVolumeFrontViewCells;
        backData.mVolume = sc.mVolumeBackViewCells;

        frontData.mCorrectedVolume = sc.mCorrectedFrontVolume;
        backData.mCorrectedVolume = sc.mCorrectedBackVolume;


    // how often was max cost ratio missed in this branch?
        frontData.mMaxCostMisses = sc.GetMaxCostMisses();
        backData.mMaxCostMisses = sc.GetMaxCostMisses();

        // set the time stamp so the order of traversal can be reconstructed
        node->SetTimeStamp(mHierarchyManager->mTimeStamp ++);
         
        // assign the objects in sorted order
        if (mUseGlobalSorting)
        {
                AssignSortedObjects(sc, frontData, backData);
        }

        // compute new stats for the view cells which see this object, 
        // e.g. new render cost decrease after the split
        UpdateViewCells(sc);

#ifdef USE_PERFTIMER  
        mNodeTimer.Exit();
#endif
        
        // return the new interior node
        return node;
}


BvhNode *BvHierarchy::Subdivide(SplitQueue &tQueue,
                                                                SubdivisionCandidate *splitCandidate,
                                                                const bool globalCriteriaMet
                                                                ,vector<SubdivisionCandidate *> &dirtyList
                                                                )
{
#ifdef USE_PERFTIMER  
        mSubdivTimer.Entry();
#endif
        
        BvhSubdivisionCandidate *sc = 
                static_cast<BvhSubdivisionCandidate *>(splitCandidate);
        BvhTraversalData &tData = sc->mParentData;

        BvhNode *currentNode = tData.mNode;

        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
        {       
                //////////////
                //-- continue subdivision

                BvhTraversalData tFrontData;
                BvhTraversalData tBackData;
                
                // create new interior node and two leaf node
                currentNode = SubdivideNode(*sc, tFrontData, tBackData);

                // decrease the weighted average cost of the subdivisoin
                mTotalCost -= sc->GetRenderCostDecrease();
                mPvsEntries += sc->GetPvsEntriesIncr();

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


                ///////////////////////////
                //-- push the new split candidates on the queue
                
                BvhSubdivisionCandidate *frontCandidate = 
                                new BvhSubdivisionCandidate(tFrontData);
                BvhSubdivisionCandidate *backCandidate = 
                                new BvhSubdivisionCandidate(tBackData);
                
                // preprocess view cells
                AssociateViewCellsWithObjects(*tData.mSampledObjects);

                EvalSubdivisionCandidate(*frontCandidate, true, false);
                EvalSubdivisionCandidate(*backCandidate, true, false);

                CollectDirtyCandidates(sc, dirtyList, true);
                
                // release preprocessed view cells
                ReleaseViewCells(*tData.mSampledObjects);
                
                // cross reference
                tFrontData.mNode->SetSubdivisionCandidate(frontCandidate); 
                tBackData.mNode->SetSubdivisionCandidate(backCandidate);

                //cout << "f: " << frontCandidate->GetPriority() << " b: " << backCandidate->GetPriority() << endl;
                tQueue.Push(frontCandidate);
                tQueue.Push(backCandidate);
        }

        /////////////////////////////////
        //-- node is a leaf => terminate traversal

        if (currentNode->IsLeaf())
        {
                /////////////////////
                //-- store additional info
                EvaluateLeafStats(tData);
        
                // this leaf is no candidate for splitting anymore
                // => detach subdivision candidate
                tData.mNode->SetSubdivisionCandidate(NULL); 
                // detach node so we don't delete it with the traversal data
                tData.mNode = NULL;
        }
        
#ifdef USE_PERFTIMER  
        mSubdivTimer.Exit();
#endif
        
        return currentNode;
}


float BvHierarchy::EvalPriority(const BvhSubdivisionCandidate &splitCandidate,
                                                                const float renderCostDecr,
                                                                const float oldRenderCost) const
{
        float priority;

        if (mIsInitialSubdivision)
        {
                priority = (float)-splitCandidate.mParentData.mDepth;
                return priority;
        }

        BvhLeaf *leaf = splitCandidate.mParentData.mNode;

        // use urface area heuristics if no view space subdivision available. 
        // For prioritized traversal we use this formula instead
        if (mHierarchyManager->GetViewSpaceSubdivisionType() == 
                HierarchyManager::NO_VIEWSPACE_SUBDIV)
        {
                priority = EvalSahCost(leaf);
        }
        else
        {
                // take render cost of node into account 
                // otherwise danger of being stuck in a local minimum!
                priority = mRenderCostDecreaseWeight          * renderCostDecr + 
                               (1.0f - mRenderCostDecreaseWeight) * oldRenderCost;
                
                if (mHierarchyManager->mConsiderMemory)
                        priority /= ((float)splitCandidate.GetPvsEntriesIncr() + mMemoryConst);
        }

        // hack: don't allow empty splits to be taken
        if (splitCandidate.mFrontObjects.empty() || splitCandidate.mBackObjects.empty())
                priority = 0;

        return priority;
}


static float AvgRayContribution(const int pvs, const int nRays)
{
        return (float)pvs / ((float)nRays + Limits::Small);
}


static float AvgRaysPerObject(const int pvs, const int nRays)
{
        return (float)nRays / ((float)pvs + Limits::Small);
}


void BvHierarchy::EvalSubdivisionCandidate(BvhSubdivisionCandidate &splitCandidate, 
                                                                                   const bool computeSplitPlane,
                                                                                   const bool preprocessViewCells)
{
#ifdef USE_PERFTIMER  
  mPlaneTimer.Entry();
#endif

        const BvhTraversalData &tData = splitCandidate.mParentData;
        BvhLeaf *leaf = tData.mNode;

#if STORE_VIEWCELLS_WITH_BVH
        if (preprocessViewCells) // fill view cells cache
                AssociateViewCellsWithObjects(*splitCandidate.mParentData.mSampledObjects);
#endif

        if (computeSplitPlane)
        {
                splitCandidate.mFrontObjects.clear();
                splitCandidate.mBackObjects.clear();

                splitCandidate.mSampledFrontObjects.clear();
                splitCandidate.mSampledBackObjects.clear();

                const bool sufficientSamples = 
                        tData.mNumRays > mMinRaysForVisibility;

                const bool useVisibiliyBasedHeuristics = 
                                        !mUseSah &&
                                        (mHierarchyManager->GetViewSpaceSubdivisionType() == 
                                        HierarchyManager::KD_BASED_VIEWSPACE_SUBDIV) &&
                                        sufficientSamples;

                // compute best object partition
                const float ratio =     SelectObjectPartition(tData, 
                                                                                                  splitCandidate.mFrontObjects, 
                                                                                                  splitCandidate.mBackObjects,
                                                                                                  useVisibiliyBasedHeuristics);
        
                // cost ratio violated?
                const bool maxCostRatioViolated = mTermMaxCostRatio < ratio;
                const int previousMisses = splitCandidate.mParentData.mMaxCostMisses;

                splitCandidate.SetMaxCostMisses(maxCostRatioViolated ? 
                                                                                previousMisses + 1 : previousMisses);

                StoreSampledObjects(splitCandidate.mSampledFrontObjects, splitCandidate.mFrontObjects);
                StoreSampledObjects(splitCandidate.mSampledBackObjects, splitCandidate.mBackObjects);
        }

#ifdef USE_PERFTIMER  
        mPlaneTimer.Exit();
#endif

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

#ifdef USE_PERFTIMER  
        mEvalTimer.Entry();
#endif
        
        // mark view cells according to what part of the split they see
        // and compute volume
        ViewCellContainer viewCells, frontViewCells, backViewCells;
        
        CollectViewCells(*tData.mSampledObjects, viewCells, false, false);
        CollectViewCells(splitCandidate.mSampledFrontObjects, frontViewCells, false, false);
        CollectViewCells(splitCandidate.mSampledBackObjects, backViewCells, false, false);

        float volFront = 0, volBack = 0, parentVol = 0;

        ViewCell::NewMail(3);

        ViewCellContainer::const_iterator fvit, fvit_end = frontViewCells.end();

        for (fvit = frontViewCells.begin(); fvit != fvit_end; ++ fvit)
        {
                ViewCell *vc = *fvit;
                vc->Mail(0);
                
                volFront += vc->GetVolume();
                parentVol += vc->GetVolume();
        }

        ViewCellContainer::const_iterator bvit, bvit_end = backViewCells.end();
        
        int frontAndBackViewCells = 0;

        for (bvit = backViewCells.begin(); bvit != bvit_end; ++ bvit)
        {
                ViewCell *vc = *bvit;

                if (vc->Mailed(0))
                {
                        // view cell sees front AND back object
                        ++ frontAndBackViewCells;
                        vc->Mail(2);
                }
                else
                {
                        vc->Mail(1);
                        parentVol += vc->GetVolume();
                }

                volBack += vc->GetVolume();
        }


        /////////////////////
        //-- this bvh node is a pvs entry in all the view cells that see one of the objects.
        
        // pvs size induced by this bvh node is #view cells
        const float pvs = (float)viewCells.size();
        
        // for low #rays per object => the result is influenced by undersampling
        const float avgRaysPerObject = AvgRaysPerObject((int)pvs, tData.mNumRays);
        splitCandidate.SetAvgRaysPerObject(avgRaysPerObject);

        const float viewSpaceVol = GetViewSpaceVolume();

        splitCandidate.mVolumeFrontViewCells = volFront / viewSpaceVol;
        splitCandidate.mVolumeBackViewCells = volBack / viewSpaceVol;

        splitCandidate.mNumFrontViewCells = (int)frontViewCells.size();
        splitCandidate.mNumBackViewCells = (int)backViewCells.size();

        
        ////////////////////////
        // warning: currently not working for new evaluation method!

        // todo matt: fix this to cope with undersampling
        splitCandidate.mCorrectedFrontVolume = 
                mHierarchyManager->EvalCorrectedPvs(splitCandidate.mVolumeFrontViewCells, 
                                                                                        parentVol, 
                                                                                        avgRaysPerObject);
        
        splitCandidate.mCorrectedBackVolume = 
                mHierarchyManager->EvalCorrectedPvs(splitCandidate.mVolumeBackViewCells, 
                                                                                        parentVol, 
                                                                                        avgRaysPerObject);

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


        float newRenderCost = 0, oldRenderCost = 0;

        // total #triangles in parent node
        const int totalTri = (int)(*tData.mSortedObjects[0]).size();
        const int frontTri = (int)splitCandidate.mFrontObjects.size();
        const int backTri = (int)splitCandidate.mBackObjects.size();
        

        // compute render cost decrease in the view cells which can see the object
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
        
        for (vit = viewCells.begin(); vit != viewCells.end(); ++ vit)
        {
                ViewCell *vc = *vit;

                const int oldVcTri = (int)vc->GetTrianglesInPvs();
                const int oldVcObj = vc->GetEntriesInPvs();

                // triangles in this view cell
                int vcTri;
                // #entries in this view cell
                int vcObj;
                
                // both nodes in this view cell
                if (vc->Mailed(2))
                {
                        vcTri = oldVcTri;
                        // #entries is increasing
                        vcObj = oldVcObj + 1;   
                }
                else if (vc->Mailed(1))
                {
                        // only back node in this view cell: #triangles is decreasing
                        vcTri = oldVcTri + backTri - totalTri; 
                        vcObj = oldVcObj;   
                }
                else // (vc->Mailed(0))
                {
                        // only front node in this view cell: #triangles is decreasing
                        vcTri = oldVcTri + frontTri - totalTri; 
                        vcObj = oldVcObj;
                }

                const float oldRc = mViewCellsManager->ComputeRenderCost(oldVcTri, oldVcObj);
                const float newRc = mViewCellsManager->ComputeRenderCost(vcTri, vcObj);

                // compute weighted render cost 
                oldRenderCost += oldRc * vc->GetVolume() / viewSpaceVol;
                newRenderCost += newRc * vc->GetVolume() / viewSpaceVol;
        }


        // compute global decrease in render cost
        const float renderCostDecr = oldRenderCost - newRenderCost;

        // for each view cell seeing both front and back object there is a new pvs entry
        splitCandidate.SetPvsEntriesIncr(frontAndBackViewCells);
        splitCandidate.SetRenderCostDecrease(renderCostDecr);
        
        const float pseudoOldRenderCost = parentVol * (float)leaf->mObjects.size() / viewSpaceVol;

        // at last computed priority based on render cost reduction / memory increase
        const float priority = EvalPriority(splitCandidate, renderCostDecr,     pseudoOldRenderCost);
        splitCandidate.SetPriority(priority);

#if STORE_VIEWCELLS_WITH_BVH
        if (preprocessViewCells)
                ReleaseViewCells(*splitCandidate.mParentData.mSampledObjects);
#endif

#ifdef USE_PERFTIMER  
        mEvalTimer.Exit();
#endif
}


int BvHierarchy::EvalPvsEntriesIncr(BvhSubdivisionCandidate &splitCandidate, 
                                                                        const float avgRaysPerObjects,
                                                                        const int numParentViewCells,
                                                                        const int numFrontViewCells,
                                                                        const int numBackViewCells) //const
{
        const float oldPvsSize = (float)numParentViewCells;
        const float oldPvsRatio = 
                (splitCandidate.mParentData.mPvs > 0) ? oldPvsSize / splitCandidate.mParentData.mPvs : 1;

        const float parentPvs = splitCandidate.mParentData.mCorrectedPvs * oldPvsRatio;

        const int frontViewCells = numFrontViewCells;
        const int backViewCells = numBackViewCells;
        
        splitCandidate.mCorrectedFrontPvs = 
                mHierarchyManager->EvalCorrectedPvs((float)frontViewCells, parentPvs, avgRaysPerObjects);
        splitCandidate.mCorrectedBackPvs =
                mHierarchyManager->EvalCorrectedPvs((float)backViewCells, parentPvs, avgRaysPerObjects);

#if GTP_DEBUG
        Debug << "bvh node pvs"
                  << " avg ray contri: " << avgRaysPerObjects << " ratio: " << oldPvsRatio 
                  << " parent: " << parentPvs << " " << " old vol: " << oldPvsSize
                  << " frontpvs: " << frontViewCells << " corr. " << splitCandidate.mCorrectedFrontPvs 
                  << " backpvs: " << frontViewCells << " corr. " << splitCandidate.mCorrectedBackPvs << endl;
#endif

        return (int)(splitCandidate.mCorrectedFrontPvs + splitCandidate.mCorrectedBackPvs - parentPvs);
}


inline bool BvHierarchy::LocalTerminationCriteriaMet(const BvhTraversalData &tData) const
{
        const bool terminationCriteriaMet =
                        (0
                        || ((int)tData.mNode->mObjects.size() <= 1)//mTermMinObjects)
                        //|| (data.mProbability <= mTermMinProbability)
                        //|| (data.mNumRays <= mTermMinRays)
                 );

#ifdef _DEBUG
        if (terminationCriteriaMet)
        {
                Debug << "bvh local termination criteria met:" << endl;
                Debug << "objects: " << (int)tData.mNode->mObjects.size() << " (" << mTermMinObjects << ")" << endl;
        }
#endif
        return terminationCriteriaMet; 
}


inline bool BvHierarchy::GlobalTerminationCriteriaMet(const BvhTraversalData &data) const
{
        // note: tracking for global cost termination 
        // does not make much sense for interleaved vsp / osp partition
        // as it is the responsibility of the hierarchy manager

        const bool terminationCriteriaMet =
                (0
                || (mBvhStats.Leaves() >= mTermMaxLeaves)
                //|| (mBvhStats.mGlobalCostMisses >= mTermGlobalCostMissTolerance)
                //|| mOutOfMemory 
                );

#ifdef GTP_DEBUG
        if (terminationCriteriaMet)
        {
                Debug << "bvh global termination criteria met:" << endl;
                Debug << "cost misses: " << mBvhStats.mGlobalCostMisses << " " << mTermGlobalCostMissTolerance << endl;
                Debug << "leaves: " << mBvhStats.Leaves() << " " << mTermMaxLeaves << endl;
        }
#endif
        return terminationCriteriaMet; 
}


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

        ////////////////
        // depth related stuff

        if (data.mDepth < mBvhStats.minDepth)
        {
                mBvhStats.minDepth = data.mDepth;
        }

        if (data.mDepth >= mTermMaxDepth)
        {
        ++ mBvhStats.maxDepthNodes;
        }

        // accumulate depth to compute average depth
        mBvhStats.accumDepth += data.mDepth;


        //////////////////////
        // objects related stuff

        // note: the sum should alwaysbe total number of objects for bvh
        mBvhStats.objectRefs += (int)leaf->mObjects.size();

        if ((int)leaf->mObjects.size() <= mTermMinObjects)
        {
             ++ mBvhStats.minObjectsNodes;
        }

        if (leaf->mObjects.empty())
        {
                ++ mBvhStats.emptyNodes;
        }

        if ((int)leaf->mObjects.size() > mBvhStats.maxObjectRefs)
        {
                mBvhStats.maxObjectRefs = (int)leaf->mObjects.size();
        }

        if ((int)leaf->mObjects.size() < mBvhStats.minObjectRefs)
        {
                mBvhStats.minObjectRefs = (int)leaf->mObjects.size();
        }

        ////////////////////////////////////////////
        // ray related stuff

        // note: this number should always accumulate to the total number of rays
        mBvhStats.rayRefs += data.mNumRays;
        
        if (data.mNumRays <= mTermMinRays)
        {
             ++ mBvhStats.minRaysNodes;
        }

        if (data.mNumRays > mBvhStats.maxRayRefs)
        {
                mBvhStats.maxRayRefs = data.mNumRays;
        }

        if (data.mNumRays < mBvhStats.minRayRefs)
        {
                mBvhStats.minRayRefs = data.mNumRays;
        }

#ifdef _DEBUG
        Debug << "depth: " << data.mDepth << " objects: " << (int)leaf->mObjects.size() 
                  << " rays: " << data.mNumRays << " rays / objects " 
                  << (float)data.mNumRays / (float)leaf->mObjects.size() << endl;
#endif
}


#if 1

/// compute object boundaries using spatial mid split
float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
                                                                                        const int axis,
                                                                                        ObjectContainer &objectsFront,
                                                                                        ObjectContainer &objectsBack)
{
        AxisAlignedBox3 parentBox = tData.mNode->GetBoundingBox();

        const float maxBox = parentBox.Max(axis);
        const float minBox = parentBox.Min(axis);

        float midPoint = (maxBox + minBox) * 0.5f;

        ObjectContainer::const_iterator oit, oit_end = tData.mNode->mObjects.end();
        
        for (oit = tData.mNode->mObjects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;
                const AxisAlignedBox3 box = obj->GetBox();

                const float objMid = (box.Max(axis) + box.Min(axis)) * 0.5f;

                // object mailed => belongs to back objects
                if (objMid < midPoint) 
                {
                        objectsBack.push_back(obj);
                }
                else
                {
                        objectsFront.push_back(obj);
                }
        }

        AxisAlignedBox3 fbox = EvalBoundingBox(objectsFront, &parentBox);
        AxisAlignedBox3 bbox = EvalBoundingBox(objectsBack, &parentBox);

        const float oldRenderCost = (float)tData.mNode->mObjects.size() * parentBox.SurfaceArea();
        const float newRenderCost = (float)objectsFront.size() * fbox.SurfaceArea() + (float)objectsBack.size() * bbox.SurfaceArea();

        const float ratio = newRenderCost / oldRenderCost;
        return ratio;
}

#else

/// compute object partition by getting balanced objects on the left and right side
float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
                                                                                        const int axis,
                                                                                        ObjectContainer &objectsFront,
                                                                                        ObjectContainer &objectsBack)
{
        PrepareLocalSubdivisionCandidates(tData, axis);
        
        SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();

        int i = 0;
        const int border = (int)tData.mNode->mObjects.size() / 2;

    for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ i)
        {
                Intersectable *obj = (*cit).mObject;

                // object mailed => belongs to back objects
                if (i < border) 
                {
                        objectsBack.push_back(obj);
                }
                else
                {
                        objectsFront.push_back(obj);
                }
        }

#if 1
        // hack: always take driving axis
        const float cost = (tData.mNode->GetBoundingBox().Size().DrivingAxis() == axis) ? -1.0f : 0.0f;
#else
        const float oldRenderCost = EvalAbsCost(tData.mLeaf->mObjects) / EvalProbability(tData.mSampledObjects);
        const float newRenderCost = EvalRenderCost(objectsFront) + EvalRenderCost(objectsBack);

        const float cost = newRenderCost / oldRenderCost;
#endif

        return cost;
}
#endif


float BvHierarchy::EvalSah(const BvhTraversalData &tData,
                                                   const int axis,
                                                   ObjectContainer &objectsFront,
                                                   ObjectContainer &objectsBack)
{
        // 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
        PrepareLocalSubdivisionCandidates(tData, axis);

        const float totalRenderCost = (float)tData.mNode->mObjects.size();
        float objectsLeft = 0, objectsRight = totalRenderCost;
  
        const AxisAlignedBox3 nodeBbox = tData.mNode->GetBoundingBox();
        const float boxArea = nodeBbox.SurfaceArea();

        float minSum = 1e20f;
  
        float minBorder = nodeBbox.Max(axis);
        float maxBorder = nodeBbox.Min(axis);

        float areaLeft = 0, areaRight = 0;

        SortableEntryContainer::const_iterator currentPos =
                mSubdivisionCandidates->begin();
        
        vector<float> bordersRight;

        // we keep track of both borders of the bounding boxes =>
        // store the events in descending order

        bordersRight.resize(mSubdivisionCandidates->size());

        SortableEntryContainer::reverse_iterator rcit = 
                mSubdivisionCandidates->rbegin(), rcit_end = mSubdivisionCandidates->rend();

        vector<float>::reverse_iterator rbit = bordersRight.rbegin();

        for (; rcit != rcit_end; ++ rcit, ++ rbit) 
        {
                Intersectable *obj = (*rcit).mObject;
                const AxisAlignedBox3 obox = obj->GetBox();

                if (obox.Min(axis) < minBorder)
                {
                        minBorder = obox.Min(axis);
                }

                (*rbit) = minBorder;
        }

        // record surface areas during the sweep
        float al = 0;
        float ar = boxArea;

        vector<float>::const_iterator bit = bordersRight.begin();
        SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();

        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit) 
        {
                Intersectable *obj = (*cit).mObject;

                ++ objectsLeft;
                -- objectsRight;

                const bool noValidSplit = ((objectsLeft <= Limits::Small) || (objectsRight <= Limits::Small));
                const AxisAlignedBox3 obox = obj->GetBox();

                // the borders of the bounding boxes have changed
                if (obox.Max(axis) > maxBorder)
                {
                        maxBorder = obox.Max(axis);
                }

                minBorder = (*bit);

                AxisAlignedBox3 lbox = nodeBbox;
                AxisAlignedBox3 rbox = nodeBbox;

                lbox.SetMax(axis, maxBorder);
                rbox.SetMin(axis, minBorder);

                al = lbox.SurfaceArea();
                ar = rbox.SurfaceArea();

                // should use classical approach here ...
#if BOUND_RENDERCOST
                const float rcLeft = std::max(objectsLeft, MIN_RENDERCOST);
                const float rcRight = std::max(objectsRight, MIN_RENDERCOST);

                const float sum = noValidSplit ? 1e25f : objectsLeft * al + objectsRight * ar;
#else

                const float sum = noValidSplit ? 1e25f : objectsLeft * al + objectsRight * ar;
#endif
        
                if (sum < minSum) 
                {       
                        minSum = sum;
                        areaLeft = al;
                        areaRight = ar;

                        // objects belong to left side now
                        for (; currentPos != (cit + 1); ++ currentPos);
                }
        }

        ////////////
        //-- assign object to front and back volume

        // belongs to back bv
        for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
                objectsBack.push_back((*cit).mObject);
        
        // belongs to front bv
        for (cit = currentPos; cit != cit_end; ++ cit)
                objectsFront.push_back((*cit).mObject);

        float newCost = minSum / boxArea;
        float ratio = newCost / totalRenderCost;
  
#ifdef GTP_DEBUG
        Debug << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of " 
                  << (int)tData.mNode->mObjects.size() << ")\t area=(" 
                  << areaLeft << ", " << areaRight << ", " << boxArea << ")" << endl
                  << "cost= " << newCost << " oldCost=" << totalRenderCost / boxArea << endl;
#endif

        return ratio;
}



float BvHierarchy::EvalSahWithTigherBbox(const BvhTraversalData &tData,
                                                                                 const int axis,
                                                                                 ObjectContainer &objectsFront,
                                                                                 ObjectContainer &objectsBack)
{
        // 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
        PrepareLocalSubdivisionCandidates(tData, axis);

        const float totalRenderCost = (float)tData.mNode->mObjects.size();
        float objectsLeft = 0, objectsRight = totalRenderCost;
  
        const AxisAlignedBox3 nodeBbox = tData.mNode->GetBoundingBox();

        const float minBox = nodeBbox.Min(axis);
        const float maxBox = nodeBbox.Max(axis);
        const float boxArea = nodeBbox.SurfaceArea();

        float minSum = 1e20f;
  
        Vector3 minBorder = nodeBbox.Max();
        Vector3 maxBorder = nodeBbox.Min();

        float areaLeft = 0, areaRight = 0;

        SortableEntryContainer::const_iterator currentPos =
                mSubdivisionCandidates->begin();
        
        vector<Vector3> bordersRight;

        // we keep track of both borders of the bounding boxes =>
        // store the events in descending order
        bordersRight.resize(mSubdivisionCandidates->size());

        SortableEntryContainer::reverse_iterator rcit = 
                mSubdivisionCandidates->rbegin(), rcit_end = 
                mSubdivisionCandidates->rend();

        vector<Vector3>::reverse_iterator rbit = bordersRight.rbegin();

        for (; rcit != rcit_end; ++ rcit, ++ rbit) 
        {
                Intersectable *obj = (*rcit).mObject;
                const AxisAlignedBox3 obox = obj->GetBox();

                for (int i = 0; i < 3; ++ i)
                {
                        if (obox.Min(i) < minBorder[i])
                        {
                                minBorder[i] = obox.Min(i);
                        }
                }

                (*rbit) = minBorder;
        }

        // temporary surface areas
        float al = 0;
        float ar = boxArea;

        vector<Vector3>::const_iterator bit = bordersRight.begin();
        SortableEntryContainer::const_iterator cit, cit_end =
                mSubdivisionCandidates->end();

        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit) 
        {
                Intersectable *obj = (*cit).mObject;

                objectsLeft ++;;
                objectsRight --;

                const AxisAlignedBox3 obox = obj->GetBox();

                AxisAlignedBox3 lbox = nodeBbox;
                AxisAlignedBox3 rbox = nodeBbox;
        
                // the borders of the left bounding box have changed
                for (int i = 0; i < 3; ++ i)
                {
                        if (obox.Max(i) > maxBorder[i])
                        {
                                maxBorder[i] = obox.Max(i);
                        }
                }

                minBorder = (*bit);

                lbox.SetMax(maxBorder);
                rbox.SetMin(minBorder);

                al = lbox.SurfaceArea();
                ar = rbox.SurfaceArea();
        
                const bool noValidSplit = ((objectsLeft <= Limits::Small) || (objectsRight <= Limits::Small));
                const float sum =  noValidSplit ? 1e25f : objectsLeft * al + objectsRight * ar;
      
                if (sum < minSum) 
                {       
                        minSum = sum;
                        areaLeft = al;
                        areaRight = ar;

                        // objects belong to left side now
                        for (; currentPos != (cit + 1); ++ currentPos);
                }
        }

        /////////////
        //-- assign object to front and back volume

        // belongs to back bv
        for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
                objectsBack.push_back((*cit).mObject);
        
        // belongs to front bv
        for (cit = currentPos; cit != cit_end; ++ cit)
                objectsFront.push_back((*cit).mObject);

        float newCost = minSum / boxArea;
        float ratio = newCost / totalRenderCost;
  
#ifdef GTP_DEBUG
        Debug << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of " 
                  << (int)tData.mNode->mObjects.size() << ")\t area=(" 
                  << areaLeft << ", " << areaRight << ", " << boxArea << ")" << endl
                  << "cost= " << newCost << " oldCost=" << totalRenderCost / boxArea << endl;
#endif

        return ratio;
}


static bool PrepareOutput(const int axis,
                                                  const int leaves,
                                                  ofstream &sumStats, 
                                                  ofstream &vollStats, 
                                                  ofstream &volrStats)
{
        if ((axis == 0) && (leaves > 0) && (leaves < 90))
        {
                char str[64];   
                sprintf(str, "tmp/bvh_heur_sum-%04d.log", leaves);
                sumStats.open(str);
                sprintf(str, "tmp/bvh_heur_voll-%04d.log", leaves);
                vollStats.open(str);
                sprintf(str, "tmp/bvh_heur_volr-%04d.log", leaves);
                volrStats.open(str);
        }

        return sumStats.is_open() && vollStats.is_open() && volrStats.is_open();
}


static void PrintHeuristics(const float objectsRight,
                                                        const float sum,
                                                        const float volLeft,
                                                        const float volRight,
                                                        const float viewSpaceVol,
                                                        ofstream &sumStats, 
                                                        ofstream &vollStats, 
                                                        ofstream &volrStats)
{
        sumStats
                << "#Position\n" << objectsRight << endl
                << "#Sum\n" << sum / viewSpaceVol << endl
                << "#Vol\n" << (volLeft +  volRight) / viewSpaceVol << endl;

        vollStats
                << "#Position\n" << objectsRight << endl
                << "#Vol\n" << volLeft / viewSpaceVol << endl;

        volrStats
                << "#Position\n" << objectsRight << endl
                << "#Vol\n" << volRight / viewSpaceVol << endl;
}


float BvHierarchy::EvalLocalCostHeuristics(const BvhTraversalData &tData,
                                                                                   const int axis,
                                                                                   ObjectContainer &objectsFront,
                                                                                   ObjectContainer &objectsBack)
{
        ////////
        // traverse split candidates, count the number of objects 
        // left and right and evaluate the cost funcion

        // prepare the heuristics, set mailboxes and counters
        const float totalVol = PrepareHeuristics(tData, axis);
        
        // local helper variables
        float volLeft = 0;
        float volRight = totalVol;
        
        const float nTotalObjects = (float)tData.mNode->mObjects.size();
        float nObjectsLeft = 0;
        float nObjectsRight = nTotalObjects;

        const float viewSpaceVol = 
                mViewCellsManager->GetViewSpaceBox().GetVolume();

        SortableEntryContainer::const_iterator backObjectsStart = 
                mSubdivisionCandidates->begin();

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

        float volBack = volLeft;
        float volFront = volRight;
        float newRenderCost = nTotalObjects * totalVol;

#ifdef GTP_DEBUG
        ofstream sumStats;
        ofstream vollStats;
        ofstream volrStats;

        const bool printStats = PrepareOutput(axis, 
                                                                                  mBvhStats.Leaves(), 
                                                                                  sumStats, 
                                                                                  vollStats, 
                                                                                  volrStats);
#endif

        ///////////////////////
        //-- the sweep heuristics
        //-- traverse through events and find best split plane

        SortableEntryContainer::const_iterator cit, 
                cit_end = cit_end = mSubdivisionCandidates->end();

        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit)
        {
                Intersectable *object = (*cit).mObject;
        
                // evaluate change in l and r volume
                // voll = view cells that see only left node (i.e., left pvs)
                // volr = view cells that see only right node (i.e., right pvs)
                EvalHeuristicsContribution(object, volLeft, volRight);

                ++ nObjectsLeft;
                -- nObjectsRight;
        
                // split is only valid if #objects on left and right is not zero
                const bool noValidSplit = (nObjectsRight <= Limits::Small);

                // the heuristics
            const float sum = noValidSplit ? 
                        1e25f : volLeft * (float)nObjectsLeft + volRight * (float)nObjectsRight;

                
#ifdef GTP_DEBUG
                if (printStats)
                {
                        PrintHeuristics(nObjectsRight, sum, volLeft, 
                                                        volRight, viewSpaceVol,
                                                        sumStats, vollStats, volrStats);
                }
#endif

                if (sum < newRenderCost)
                {
                        newRenderCost = sum;

                        volBack = volLeft;
                        volFront = volRight;

                        // objects belongs to left side now
                        for (; backObjectsStart != (cit + 1); ++ backObjectsStart);
                }
        }

        ////////////////////////////////////////
        //-- assign object to front and back volume

        // belongs to back bv
        for (cit = mSubdivisionCandidates->begin(); cit != backObjectsStart; ++ cit)
        {
                objectsBack.push_back((*cit).mObject);
        }
        // belongs to front bv
        for (cit = backObjectsStart; cit != cit_end; ++ cit)
        {
                objectsFront.push_back((*cit).mObject);
        }

        // render cost of the old parent
        const float oldRenderCost = (float)nTotalObjects * totalVol + Limits::Small;
        // the relative cost ratio
        const float ratio = newRenderCost / oldRenderCost;

#ifdef GTP_DEBUG
        Debug << "\neval bvh split cost decrease" << endl
                  << "back pvs: " << (int)objectsBack.size() << " front pvs: " 
                  << (int)objectsFront.size() << " total pvs: " << nTotalObjects << 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;
#endif

        return ratio;
}


void BvHierarchy::PrepareLocalSubdivisionCandidates(const BvhTraversalData &tData,
                                                                                                        const int axis)                                                                                 
{
#ifdef USE_PERFTIMER  
        mSortTimer.Entry();
#endif  
        //-- insert object queries
        ObjectContainer *objects = mUseGlobalSorting ? 
                tData.mSortedObjects[axis] : &tData.mNode->mObjects;

        CreateLocalSubdivisionCandidates(*objects, &mSubdivisionCandidates, !mUseGlobalSorting, axis);
        
#ifdef USE_PERFTIMER  
        mSortTimer.Exit();
#endif
}


void BvHierarchy::CreateLocalSubdivisionCandidates(const ObjectContainer &objects, 
                                                                                                  SortableEntryContainer **subdivisionCandidates, 
                                                                                                  const bool sortEntries,
                                                                                                  const int axis)
{
        (*subdivisionCandidates)->clear();

        // compute requested size and look if subdivision candidate has to be recomputed
        const int requestedSize = (int)objects.size();
        
        // creates a sorted split candidates array
        if ((*subdivisionCandidates)->capacity() > 500000 &&
                requestedSize < (int)((*subdivisionCandidates)->capacity() / 10) )
        {
        delete (*subdivisionCandidates);
                (*subdivisionCandidates) = new SortableEntryContainer;
        }

        (*subdivisionCandidates)->reserve(requestedSize);

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

        for (oit = objects.begin(); oit < oit_end; ++ oit)
        {
                (*subdivisionCandidates)->push_back(SortableEntry(*oit, (*oit)->GetBox().Center(axis)));
        }

        if (sortEntries)
        {       // no presorted candidate list
                stable_sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
                //sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
        }
}


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


float BvHierarchy::PrepareHeuristics(const BvhTraversalData &tData, 
                                                                         const int axis)
{       
        BvhLeaf *leaf = tData.mNode;
        float vol = 0;

    // sort so we can use a sweep from right to left
        PrepareLocalSubdivisionCandidates(tData, axis);
        
        // collect and mark the view cells as belonging to front pvs
        ViewCellContainer viewCells;

        const bool setCounter = true;
        const bool onlyUnmailed = true;

        
        CollectViewCells(*tData.mSampledObjects, 
                                         viewCells, 
                                         setCounter, 
                                         onlyUnmailed);

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

        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
        {
#if USE_VOLUMES_FOR_HEURISTICS
                const float volIncr = (*vit)->GetVolume();
#else
                const float volIncr = 1.0f;
#endif
                vol += volIncr;
        }

        // mail view cells that go from front node to back node
        ViewCell::NewMail();
        
        return vol;
}



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


void BvHierarchy::EvalHeuristicsContribution(Intersectable *obj, 
                                                                                         float &volLeft, 
                                                                                         float &volRight)
{
        // collect all view cells associated with this objects 
        // (also multiple times, if they are pierced by several rays)
        ViewCellContainer viewCells;

        const bool useMailboxing = false;
        const bool setCounter = false;
        const bool onlyUnmailedRays = true;

        CollectViewCells(obj, viewCells, useMailboxing, setCounter, onlyUnmailedRays);

        // 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;

#if USE_VOLUMES_FOR_HEURISTICS
                const float vol = viewCell->GetVolume();
#else
                const float vol = 1.0f;
#endif
                if (!viewCell->Mailed())
                {
                        viewCell->Mail();
                        // we now see view cell from both nodes 
                        // => add volume to left node
                        volLeft += vol;
                }

                // last reference into the right node
                if (-- viewCell->mCounter == 0)
                {       
                        // view cell was previously seen from both nodes  =>
                        // remove volume from right node
                        volRight -= vol;
                }
        }
}


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


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


float BvHierarchy::SelectObjectPartition(const BvhTraversalData &tData,
                                                                                 ObjectContainer &frontObjects,
                                                                                 ObjectContainer &backObjects,
                                                                                 bool useVisibilityBasedHeuristics)
{
#ifdef USE_PERFTIMER  
        mSplitTimer.Entry();
#endif
        
        if (mIsInitialSubdivision)
        {
                ApplyInitialSplit(tData, frontObjects, backObjects);
                return 0;
        }

        ObjectContainer nFrontObjects[3];
        ObjectContainer nBackObjects[3];
        float nCostRatio[3];

        int sAxis = 0;
        int bestAxis = -1;

        if (mOnlyDrivingAxis)
        {
                const AxisAlignedBox3 box = tData.mNode->GetBoundingBox();
                sAxis = box.Size().DrivingAxis();
        }

        // if #rays high, consider only use a subset of the rays for 
        // visibility based heuristics
        VssRay::NewMail();


        ////////////////////////////////////
        //-- evaluate split cost for all three axis
        
        for (int axis = 0; axis < 3; ++ axis)
        {
                if (!mOnlyDrivingAxis || (axis == sAxis))
                {
                        if (mUseCostHeuristics)
                        {
                                //////////////////////////////////
                //-- split objects using heuristics
                                
                                if (useVisibilityBasedHeuristics)
                                {
                                        ///////////
                                        //-- heuristics using objects weighted by view cells volume
                                        nCostRatio[axis] =
                                                EvalLocalCostHeuristics(tData, 
                                                                                                axis, 
                                                                                                nFrontObjects[axis], 
                                                                                                nBackObjects[axis]);
                                }
                                else
                                {       
                                        //////////////////
                                        //-- view cells not constructed yet     => use surface area heuristic                   
                                        nCostRatio[axis] = EvalSah(tData, 
                                                                                           axis, 
                                                                                           nFrontObjects[axis], 
                                                                                           nBackObjects[axis]);
                                }
                        }
                        else
                        {
                                //-- split objects using some simple criteria
                                nCostRatio[axis] =
                                        EvalLocalObjectPartition(tData, axis, nFrontObjects[axis], nBackObjects[axis]);
                        }

                        // avoid splits in degenerate axis with high penalty
                        if (1 &&
                                (tData.mNode->GetBoundingBox().Size(axis) < 0.0001))//Limits::Small))
                        {
                                nCostRatio[axis] += 9999;
                        }

                        if ((bestAxis == -1) || (nCostRatio[axis] < nCostRatio[bestAxis]))
                        {
                                bestAxis = axis;
                        }
                }
        }

    ////////////////
        //-- assign values

        frontObjects = nFrontObjects[bestAxis];
        backObjects = nBackObjects[bestAxis];

#ifdef USE_PERFTIMER  
        mSplitTimer.Exit();
#endif
        
        //cout << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
        return nCostRatio[bestAxis];
}


int BvHierarchy::AssociateObjectsWithRays(const VssRayContainer &rays) const
{
        int nRays = 0;
        VssRayContainer::const_iterator rit, rit_end = rays.end();

        VssRay *lastVssRay = NULL;

        VssRay::NewMail();

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

                // filter out double rays (last ray the same as this ray)
                if (
                        !lastVssRay ||
                        !(ray->mOrigin == lastVssRay->mTermination) ||
                        !(ray->mTermination == lastVssRay->mOrigin)) 
                {
                        lastVssRay = ray;
                        //cout << "j";
                        if (ray->mTerminationObject)
                        {
                                ray->mTerminationObject->GetOrCreateRays()->push_back(ray);
                                if (!ray->Mailed())
                                {
                                        ray->Mail();
                                        ++ nRays;
                                }
                        }

#if COUNT_ORIGIN_OBJECTS

                        if (ray->mOriginObject)
                        {
                                //cout << "o";
                                ray->mOriginObject->GetOrCreateRays()->push_back(ray);

                                if (!ray->Mailed())
                                {
                                        ray->Mail();
                                        ++ nRays;
                                }
                        }
#endif
                }
        }

        return nRays;
}


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

        mSubdivisionStats 
                        << "#Leaves\n" << mBvhStats.Leaves() << endl
                        << "#RenderCostDecrease\n" << costDecr << endl 
                        << "#TotalRenderCost\n" << mTotalCost << endl
                        << "#EntriesInPvs\n" << mPvsEntries << 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->GetOrCreateRays()->end();

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

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


float BvHierarchy::EvalSahCost(BvhLeaf *leaf) const
{
        ////////////////
        //-- surface area heuristics

        const AxisAlignedBox3 box = GetBoundingBox(leaf);
        const float area = box.SurfaceArea();
        const float viewSpaceArea = mViewCellsManager->GetViewSpaceBox().SurfaceArea();

        return (float)leaf->mObjects.size() * area / viewSpaceArea;
}


float BvHierarchy::EvalRenderCost(const ObjectContainer &objects)// const
{       
        ///////////////
        //-- render cost heuristics

        const float objRenderCost = (float)objects.size();

        const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();

        // probability that view point lies in a view cell which sees this node
        const float p = EvalViewCellsVolume(objects) / viewSpaceVol;
        
        return objRenderCost * p;
}


float BvHierarchy::EvalProbability(const ObjectContainer &objects)// const
{       
        const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
        
        // probability that view point lies in a view cell which sees this node
        return EvalViewCellsVolume(objects) / viewSpaceVol;
}


AxisAlignedBox3 BvHierarchy::EvalBoundingBox(const ObjectContainer &objects,
                                                                                         const AxisAlignedBox3 *parentBox) const
{
        // if there are no objects in this box, box size is set to parent box size. 
        // Question: Invalidate box instead?
        if (parentBox && objects.empty())
                return *parentBox;

        AxisAlignedBox3 box;
        box.Initialize();

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

        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;
                // grow bounding box to include all objects
                box.Include(obj->GetBox());
        }

        return box;
}


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

        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());
                }
        }
}


void BvHierarchy::CollectNodes(BvhNode *root, vector<BvhNode *> &nodes) const
{
        stack<BvhNode *> nodeStack;
        nodeStack.push(root);

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

                nodes.push_back(node);
                
                if (!node->IsLeaf()) 
                {
                        BvhInterior *interior = (BvhInterior *)node;

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


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


int BvHierarchy::CollectViewCells(const ObjectContainer &objects, 
                                                                  ViewCellContainer &viewCells,
                                                                  const bool setCounter,
                                                                  const bool onlyUnmailedRays)// const
{
        ViewCell::NewMail();

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

        // use mailing to avoid dublicates
        const bool useMailBoxing = true;

        int numRays = 0;
        // loop through all object and collect view cell pvs of this node
        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                // use mailing to avoid duplicates
                numRays += CollectViewCells(*oit, viewCells, useMailBoxing, setCounter, onlyUnmailedRays);
        }

        return numRays;
}


#if STORE_VIEWCELLS_WITH_BVH


void BvHierarchy::ReleaseViewCells(const ObjectContainer &objects)
{
        ObjectContainer::const_iterator oit, oit_end = objects.end();

        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                (*oit)->DelViewCells();
        }
}


void BvHierarchy::AssociateViewCellsWithObjects(const ObjectContainer &objects) const
{
        ObjectContainer::const_iterator oit, oit_end = objects.end();

        const bool useMailBoxing = true;
        VssRay::NewMail();
        
        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                        ViewCell::NewMail();
                        // use mailing to avoid duplicates
                        AssociateViewCellsWithObject(*oit, useMailBoxing);
        }
}


int BvHierarchy::AssociateViewCellsWithObject(Intersectable *obj, const bool useMailBoxing) const
{
        int nRays = 0;

        if (!obj->GetOrCreateViewCells()->empty())
        {
                cerr << "AssociateViewCellsWithObject: view cells cache not working" << endl;
        }

        ViewCellContainer *objViewCells = obj->GetOrCreateViewCells();
        VssRayContainer *vssRays = obj->GetOrCreateRays();

        VssRayContainer::const_iterator rit, rit_end = vssRays->end();

        // fill cache
        for (rit = vssRays->begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit);

                //      if (onlyUnmailedRays && ray->Mailed())
                //              continue;
                mHierarchyManager->mVspTree->GetViewCells(*ray, *objViewCells);
                
                if (!useMailBoxing || !ray->Mailed())
                {
                        if (useMailBoxing)
                                ray->Mail();

                        ++ nRays;
                }
        }

        return nRays;
}



int BvHierarchy::CountViewCells(Intersectable *obj) //const
{
        ViewCellContainer *viewCells = obj->GetOrCreateViewCells();

        if (obj->GetOrCreateViewCells()->empty())
        {
                //cerr << "h";//CountViewCells: view cells empty, view cells cache not working" << endl;
                return CountViewCellsFromRays(obj);
        }
        
        int result = 0;

        ViewCellContainer::const_iterator vit, vit_end = viewCells->end();
        
        for (vit = viewCells->begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = *vit;

                // store view cells
                if (!vc->Mailed())
                {
                        vc->Mail();
                        ++ result;
                }
        }

        return result;
}


int BvHierarchy::CollectViewCells(Intersectable *obj, 
                                                                  ViewCellContainer &viewCells,
                                                                  const bool useMailBoxing,
                                                                  const bool setCounter,
                                                                  const bool onlyUnmailedRays)// const
{
        // view cells not cached
        if (obj->GetOrCreateViewCells()->empty())
        {
                return CollectViewCellsFromRays(obj, viewCells, useMailBoxing, setCounter, onlyUnmailedRays);
        }

        ///////////
        //-- use view cells cache

#ifdef USE_PERFTIMER  
        mCollectTimer.Entry();
#endif
        
        ViewCellContainer *objViewCells = obj->GetOrCreateViewCells();

        // loop through view cells
        // matt: probably slow to insert view cells one by one
        ViewCellContainer::const_iterator vit, vit_end = objViewCells->end();

        for (vit = objViewCells->begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = *vit;

                // store view cells
                if (!useMailBoxing || !vc->Mailed())
                {
                        if (useMailBoxing) 
                        {
                                // view cell not mailed
                                vc->Mail();
                                
                                if (setCounter)
                                        vc->mCounter = 0;
                                //viewCells.push_back(vc);
                        }

                        viewCells.push_back(vc);
                }

                if (setCounter)
                        ++ vc->mCounter;
        }

#ifdef USE_PERFTIMER  
        mCollectTimer.Exit();
#endif
        return (int)objViewCells->size();
}


int BvHierarchy::CollectViewCellsFromRays(Intersectable *obj, 
                                                                                  ViewCellContainer &viewCells,
                                                                                  const bool useMailBoxing,
                                                                                  const bool setCounter,
                                                                                  const bool onlyUnmailedRays)
{
#ifdef USE_PERFTIMER  
        mCollectTimer.Entry();
#endif
        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

        int numRays = 0;

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

                if (onlyUnmailedRays && ray->Mailed())
                        continue;
                
                ++ numRays;

                ViewCellContainer tmpViewCells;
                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);

                // matt: probably slow to allocate memory for view cells every time
                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();

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

                        // store view cells
                        if (!useMailBoxing || !vc->Mailed())
                        {
                                if (useMailBoxing) // => view cell not mailed
                                {
                                        vc->Mail();
                                        if (setCounter)
                                                vc->mCounter = 0;
                                }

                                viewCells.push_back(vc);
                        }
                        
                        if (setCounter)
                                ++ vc->mCounter;
                }
        }

#ifdef USE_PERFTIMER  
        mCollectTimer.Exit();
#endif
        return numRays;
}


int BvHierarchy::CountViewCellsFromRays(Intersectable *obj) //const
{
        int result = 0;
        
        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

        for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit);
                ViewCellContainer tmpViewCells;
        
                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
                
                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;

                        // store view cells
                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                ++ result;
                        }
                }
        }

        return result;
}

#else

int BvHierarchy::CountViewCells(Intersectable *obj) //const
{
        int result = 0;
        
        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

        for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit);
                ViewCellContainer tmpViewCells;
        
                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
                
                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;

                        // store view cells
                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                ++ result;
                        }
                }
        }

        return result;
}


int BvHierarchy::CollectViewCells(Intersectable *obj, 
                                                                  ViewCellContainer &viewCells,
                                                                  const bool useMailBoxing,
                                                                  const bool setCounter,
                                                                  const bool onlyUnmailedRays)
{
#ifdef USE_PERFTIMER  
        mCollectTimer.Entry();
#endif
        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();

        int numRays = 0;

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

                if (onlyUnmailedRays && ray->Mailed())
                        continue;
                
                ++ numRays;

                ViewCellContainer tmpViewCells;
                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);

                // matt: probably slow to allocate memory for view cells every time
                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();

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

                        // store view cells
                        if (!useMailBoxing || !vc->Mailed())
                        {
                                if (useMailBoxing) // => view cell not mailed
                                {
                                        vc->Mail();
                                        if (setCounter)
                                                vc->mCounter = 0;
                                }

                                viewCells.push_back(vc);
                        }
                        
                        if (setCounter)
                                ++ vc->mCounter;
                }
        }

#ifdef USE_PERFTIMER  
        mCollectTimer.Exit();
#endif
        return numRays;
}
#endif


int BvHierarchy::CountViewCells(const ObjectContainer &objects)// const
{
        int nViewCells = 0;

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

        // loop through all object and collect view cell pvs of this node
        for (oit = objects.begin(); oit != oit_end; ++ oit)
        {
                nViewCells += CountViewCells(*oit);
        }

        return nViewCells;
}


void BvHierarchy::CollectDirtyCandidates(BvhSubdivisionCandidate *sc, 
                                                                                 vector<SubdivisionCandidate *> &dirtyList, 
                                                                                 const bool onlyUnmailed)
{
        BvhTraversalData &tData = sc->mParentData;
        BvhLeaf *node = tData.mNode;
        
        ViewCellContainer viewCells;
        //ViewCell::NewMail();
        int numRays = CollectViewCells(*tData.mSampledObjects, viewCells, false, false);

        if (0) cout << "collected " << (int)viewCells.size() << " dirty candidates" << endl;
        
        // split candidates handling 
        // these 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 = static_cast<VspViewCell *>(*vit);
                VspLeaf *leaf = vc->mLeaves[0];
        
                SubdivisionCandidate *candidate = leaf->GetSubdivisionCandidate();
                
                // is this leaf still a split candidate?
                if (candidate && (!onlyUnmailed || !candidate->Mailed()))
                {
                        candidate->Mail();
                        candidate->SetDirty(true);
                        dirtyList.push_back(candidate);
                }
        }
}


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;
        }
}

#if 0
BvhLeaf *BvHierarchy::GetLeaf(Intersectable *object, BvhNode *node) const
{
        // hack: we use the simpler but faster version
        if (!object) 
                return NULL;

        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 = static_cast<BvhLeaf *>(node);

                        if (IsObjectInLeaf(leaf, object))
                        {
                                return leaf;
                        }
                } 
                else    
                {       
                        // find point
                        BvhInterior *interior = static_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;
}
#endif

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

        return true;
}


void BvHierarchy::ExportObjects(BvhLeaf *leaf, OUT_STREAM &stream)
{
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                stream << (*oit)->GetId() << " ";
        }
}


void BvHierarchy::ExportNode(BvhNode *node, OUT_STREAM &stream)
{
        if (node->IsLeaf())
        {
                BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
                const AxisAlignedBox3 box = leaf->GetBoundingBox();
                stream << "<Leaf id=\"" << node->GetId() << "\""
                           << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
                           << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\""
                           << " objects=\"";
                
                //-- export objects
                // tmp matt
                if (1) ExportObjects(leaf, stream);
                
                stream << "\" />" << endl;
        }
        else
        {       
                BvhInterior *interior = static_cast<BvhInterior *>(node);
                const AxisAlignedBox3 box = interior->GetBoundingBox();

                stream << "<Interior id=\"" << node->GetId() << "\""
                           << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
                           << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z
                           << "\">" << endl;

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

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


float BvHierarchy::EvalViewCellsVolume(const ObjectContainer &objects)// const
{
        float vol = 0;

        ViewCellContainer viewCells;
        
        // we have to account for all view cells that can 
        // be seen from the objects
        int numRays = CollectViewCells(objects, viewCells, false, false);

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

        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
        {
                vol += (*vit)->GetVolume();
        }

        return vol;
}


void BvHierarchy::Initialise(const ObjectContainer &objects)
{
        AxisAlignedBox3 box = EvalBoundingBox(objects);

        ///////
        //-- create new root

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

        // compute bounding box from objects
        mBoundingBox = mRoot->GetBoundingBox();

        // associate root with current objects
        AssociateObjectsWithLeaf(bvhleaf);
}


void BvHierarchy::StoreSampledObjects(ObjectContainer &sampledObjects, const ObjectContainer &objects)
{
        ObjectContainer::const_iterator oit, oit_end = objects.end();

        for (oit = objects.begin(); oit != objects.end(); ++ oit)
        {
                Intersectable *obj = *oit;

                if (!obj->GetOrCreateRays()->empty())
        {
                        sampledObjects.push_back(obj);
                }
        }
        }


void BvHierarchy::PrepareConstruction(SplitQueue &tQueue,
                                                                          const VssRayContainer &sampleRays,
                                                                          const ObjectContainer &objects)
{
        ///////////////////////////////////////
        //-- we assume that we have objects sorted by their id => 
        //-- we don't have to sort them here and an binary search 
        //-- for identifying if a object is in a leaf.
        
        mBvhStats.Reset();
        mBvhStats.Start();
        mBvhStats.nodes = 1;
                
        // store pointer to this tree
        BvhSubdivisionCandidate::sBvHierarchy = this;
        
        // root and bounding box was already constructed
        BvhLeaf *bvhLeaf = static_cast<BvhLeaf *>(mRoot);
        
        // only rays intersecting objects in node are interesting
        const int nRays = AssociateObjectsWithRays(sampleRays);
        //cout << "using " << nRays << " of " << (int)sampleRays.size() << " rays" << endl;
        
        ObjectContainer *sampledObjects = new ObjectContainer();
        StoreSampledObjects(*sampledObjects, objects);

#if STORE_VIEWCELLS_WITH_BVH
        AssociateViewCellsWithObjects(*sampledObjects);
#endif

        // probability that volume is "seen" from the view cells
        const float prop = EvalViewCellsVolume(*sampledObjects) / GetViewSpaceVolume();

        // create bvh traversal data
        BvhTraversalData oData(bvhLeaf, 0, prop, nRays);
                
        // create sorted object lists for the first data
        if (mUseGlobalSorting)
        {
                AssignInitialSortedObjectList(oData, objects);
        }
        
        oData.mSampledObjects = sampledObjects;
        
        ///////////////////
        //-- add first candidate for object space partition     

        mTotalCost = EvalRenderCost(objects);
        mPvsEntries = CountViewCells(*sampledObjects);

        oData.mCorrectedPvs = oData.mPvs = (float)mPvsEntries;
        oData.mCorrectedVolume = oData.mVolume = prop;
        
        BvhSubdivisionCandidate *oSubdivisionCandidate = 
                new BvhSubdivisionCandidate(oData);

        bvhLeaf->SetSubdivisionCandidate(oSubdivisionCandidate);

#if STORE_VIEWCELLS_WITH_BVH
        ReleaseViewCells(*sampledObjects);
#endif

        if (mApplyInitialPartition)
        {
                vector<SubdivisionCandidate *> candidateContainer;

                mIsInitialSubdivision = true;
                
                // evaluate priority
                EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
                PrintSubdivisionStats(*oSubdivisionCandidate);

                ApplyInitialSubdivision(oSubdivisionCandidate, candidateContainer);             

                mIsInitialSubdivision = false;

                vector<SubdivisionCandidate *>::const_iterator cit, cit_end = candidateContainer.end();

                for (cit = candidateContainer.begin(); cit != cit_end; ++ cit)
                {
                        BvhSubdivisionCandidate *sCandidate = static_cast<BvhSubdivisionCandidate *>(*cit);
                        
                        // reevaluate priority
                        EvalSubdivisionCandidate(*sCandidate, true, true);
                        tQueue.Push(sCandidate);
                }

                cout << "size of initial bv subdivision: " << GetStatistics().Leaves() << endl;
        }
        else
        {       
                // evaluate priority
                EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
                PrintSubdivisionStats(*oSubdivisionCandidate);

                tQueue.Push(oSubdivisionCandidate);
                cout << "size of initial bv subdivision: " << GetStatistics().Leaves() << endl;
        }
}


void BvHierarchy::AssignInitialSortedObjectList(BvhTraversalData &tData,
                                                                                                const ObjectContainer &objects)
{
        const bool doSort = true;

        // we sort the objects as a preprocess so they don't have
        // to be sorted for each split
        for (int i = 0; i < 3; ++ i)
        {
                SortableEntryContainer *sortedObjects = new SortableEntryContainer();

                CreateLocalSubdivisionCandidates(objects, 
                                                                             &sortedObjects, 
                                                                                 doSort, 
                                                                                 i);
                
                // copy list into traversal data list
                tData.mSortedObjects[i] = new ObjectContainer();
                tData.mSortedObjects[i]->reserve((int)objects.size());

                SortableEntryContainer::const_iterator oit, oit_end = sortedObjects->end();

                for (oit = sortedObjects->begin(); oit != oit_end; ++ oit)
                {
                        tData.mSortedObjects[i]->push_back((*oit).mObject);
                }

                delete sortedObjects;
        }

        // next sorted list: by size (for initial criteria)
        tData.mSortedObjects[3] = new ObjectContainer();
        tData.mSortedObjects[3]->reserve((int)objects.size());

        *(tData.mSortedObjects[3]) = objects;
        
        stable_sort(tData.mSortedObjects[3]->begin(), tData.mSortedObjects[3]->end(), smallerSize);
}


void BvHierarchy::AssignSortedObjects(const BvhSubdivisionCandidate &sc,
                                                                          BvhTraversalData &frontData,
                                                                          BvhTraversalData &backData)
{
        Intersectable::NewMail();

        // we sorted the objects as a preprocess so they don't have
        // to be sorted for each split
        ObjectContainer::const_iterator fit, fit_end = sc.mFrontObjects.end();

        for (fit = sc.mFrontObjects.begin(); fit != fit_end; ++ fit)
        {
                (*fit)->Mail();
        }

        for (int i = 0; i < 4; ++ i)
        {
                frontData.mSortedObjects[i] = new ObjectContainer();
                backData.mSortedObjects[i] = new ObjectContainer();

                frontData.mSortedObjects[i]->reserve(sc.mFrontObjects.size());
                backData.mSortedObjects[i]->reserve(sc.mBackObjects.size());

                ObjectContainer::const_iterator oit, oit_end = sc.mParentData.mSortedObjects[i]->end();

                // all the front objects are mailed => assign the sorted object lists
                for (oit = sc.mParentData.mSortedObjects[i]->begin(); oit != oit_end; ++ oit)
                {
                        if ((*oit)->Mailed())
                        {
                                frontData.mSortedObjects[i]->push_back(*oit);
                        }
                        else
                        {
                                backData.mSortedObjects[i]->push_back(*oit);
                        }
                }
        }
}


void BvHierarchy::Reset(SplitQueue &tQueue,
                                                const VssRayContainer &sampleRays, 
                                                const ObjectContainer &objects)
{

        // reset stats
        mBvhStats.Reset();
        mBvhStats.Start();
        mBvhStats.nodes = 1;

        // reset root
        DEL_PTR(mRoot);
        
        BvhLeaf *bvhleaf = new BvhLeaf(mBoundingBox, NULL, (int)objects.size());
        bvhleaf->mObjects = objects;
        mRoot = bvhleaf;
        
        ObjectContainer *sampledObjects = new ObjectContainer();
        StoreSampledObjects(*sampledObjects, objects);

#if STORE_VIEWCELLS_WITH_BVH
        AssociateViewCellsWithObjects(*sampledObjects);
#endif

        //mTermMinProbability *= mVspTree->GetBoundingBox().GetVolume();
        // probability that volume is "seen" from the view cells
        const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
        const float prop = EvalViewCellsVolume(*sampledObjects);

        const int nRays = CountRays(*sampledObjects);
        BvhLeaf *bvhLeaf = static_cast<BvhLeaf *>(mRoot);

        // create bvh traversal data
        BvhTraversalData oData(bvhLeaf, 0, prop, nRays);

        oData.mSampledObjects = sampledObjects;

        if (mUseGlobalSorting)
                AssignInitialSortedObjectList(oData, objects);
        
#if STORE_VIEWCELLS_WITH_BVH
        ReleaseViewCells(*sampledObjects);
#endif
        ///////////////////
        //-- add first candidate for object space partition     

        BvhSubdivisionCandidate *oSubdivisionCandidate = 
                new BvhSubdivisionCandidate(oData);

        EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
        bvhLeaf->SetSubdivisionCandidate(oSubdivisionCandidate);

        mTotalCost = (float)objects.size() * prop;

        PrintSubdivisionStats(*oSubdivisionCandidate);

        tQueue.Push(oSubdivisionCandidate);
}


void BvhStatistics::Print(ostream &app) const
{
        app << "=========== BvHierarchy statistics ===============\n";

        app << setprecision(4);

        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";

        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";

        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";

        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";

        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits << endl;

        app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
                << maxCostNodes * 100 / (double)Leaves() << endl;

        app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
                << minProbabilityNodes * 100 / (double)Leaves() << endl;


        //////////////////////////////////////////////////
        
        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n" 
                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
        
        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;

        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;

        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;

        
        ////////////////////////////////////////////////////////
        
        app << "#N_PMINOBJECTSLEAVES  ( Percentage of leaves with mininum objects )\n"
                << minObjectsNodes * 100 / (double)Leaves() << endl;

        app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" << maxObjectRefs << "\n";

        app << "#N_MINOBJECTREFS  ( Min number of object refs / leaf )\n" << minObjectRefs << "\n";

        app << "#N_EMPTYLEAFS ( Empty leafs )\n" << emptyNodes << "\n";
        
        app << "#N_PAVGOBJECTSLEAVES  ( average object refs / leaf)\n" << AvgObjectRefs() << endl;


        ////////////////////////////////////////////////////////
        
        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with mininum rays )\n"
                << minRaysNodes * 100 / (double)Leaves() << endl;

        app << "#N_MAXRAYREFS  ( Max number of ray refs / leaf )\n" << maxRayRefs << "\n";

        app << "#N_MINRAYREFS  ( Min number of ray refs / leaf )\n" << minRayRefs << "\n";
        
        app << "#N_PAVGRAYLEAVES  ( average ray refs / leaf )\n" << AvgRayRefs() << endl;
        
        app << "#N_PAVGRAYCONTRIBLEAVES  ( Average ray contribution)\n" <<
                rayRefs / (double)objectRefs << endl;

        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
                maxRayContriNodes * 100 / (double)Leaves() << endl;

        app << "#N_PGLOBALCOSTMISSES ( Global cost misses )\n" << mGlobalCostMisses << endl;

        app << "========== END OF BvHierarchy statistics ==========\n";
}


// TODO: return memory usage in MB
float BvHierarchy::GetMemUsage() const
{
        return (float)(sizeof(BvHierarchy)
                                   + mBvhStats.Leaves() * sizeof(BvhLeaf) 
                                   + mBvhStats.Interior() * sizeof(BvhInterior)
                                   ) / float(1024 * 1024);
}


void BvHierarchy::SetActive(BvhNode *node) const
{
        vector<BvhLeaf *> leaves;

        // sets the pointers to the currently active view cells
        CollectLeaves(node, leaves);
        vector<BvhLeaf *>::const_iterator lit, lit_end = leaves.end();

        for (lit = leaves.begin(); lit != lit_end; ++ lit)
        {
                (*lit)->SetActiveNode(node);
        }
}


void BvHierarchy::CollectObjects(const AxisAlignedBox3 &box,
                                                                 ObjectContainer &objects)
{
  stack<BvhNode *> nodeStack;
  
  nodeStack.push(mRoot);

  while (!nodeStack.empty()) {
        BvhNode *node = nodeStack.top();
        
        nodeStack.pop();
        if (node->IsLeaf()) {
          BvhLeaf *leaf = (BvhLeaf *)node;
          if (Overlap(box, leaf->GetBoundingBox())) {
                Intersectable *object = leaf;
                if (!object->Mailed()) {
                  object->Mail();
                  objects.push_back(object);
                }
          }
        } 
        else 
          {
                BvhInterior *interior = (BvhInterior *)node;
                if (Overlap(box, interior->GetBoundingBox())) {
                  bool pushed = false;
                  if (!interior->GetFront()->Mailed()) {
                        nodeStack.push(interior->GetFront());
                        pushed = true;
                  }
                  if (!interior->GetBack()->Mailed()) {
                        nodeStack.push(interior->GetBack());
                        pushed = true;
                  }
                  // avoid traversal of this node in the next query
                  if (!pushed)
                        interior->Mail();
                }
          }
  }
}


void BvHierarchy::CreateUniqueObjectIds()
{
        stack<BvhNode *> nodeStack;
        nodeStack.push(mRoot);

        int currentId = 0;
        while (!nodeStack.empty()) 
        {
                BvhNode *node = nodeStack.top();
                nodeStack.pop();

                node->SetId(currentId ++);

                if (!node->IsLeaf()) 
                {
                        BvhInterior *interior = (BvhInterior *)node;

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


void BvHierarchy::ApplyInitialSubdivision(SubdivisionCandidate *firstCandidate,
                                                                                  vector<SubdivisionCandidate *> &candidateContainer)
{
        SplitQueue tempQueue;
        tempQueue.Push(firstCandidate);

        while (!tempQueue.Empty())
        {
                SubdivisionCandidate *candidate = tempQueue.Top();
                tempQueue.Pop();

                BvhSubdivisionCandidate *bsc = 
                        static_cast<BvhSubdivisionCandidate *>(candidate);

                if (!InitialTerminationCriteriaMet(bsc->mParentData))
                {
                        const bool globalCriteriaMet = GlobalTerminationCriteriaMet(bsc->mParentData);
                
                        SubdivisionCandidateContainer dirtyList;
                        BvhNode *node = Subdivide(tempQueue, bsc, globalCriteriaMet, dirtyList);

                        // not needed anymore
                        delete bsc;
                }
                else 
                {
                        // initial preprocessing  finished for this candidate
                        // add to candidate container
                        candidateContainer.push_back(bsc);
                }
        }
}


void BvHierarchy::ApplyInitialSplit(const BvhTraversalData &tData,
                                                                        ObjectContainer &frontObjects,
                                                                        ObjectContainer &backObjects)
{
        ObjectContainer *objects = tData.mSortedObjects[3];

        ObjectContainer::const_iterator oit, oit_end = objects->end();
   
        float maxAreaDiff = -1.0f;

        ObjectContainer::const_iterator backObjectsStart = objects->begin();

        for (oit = objects->begin(); oit != (objects->end() - 1); ++ oit)
        {
                Intersectable *objS = *oit;
                Intersectable *objL = *(oit + 1);
                
                const float areaDiff = 
                                objL->GetBox().SurfaceArea() - objS->GetBox().SurfaceArea();

                if (areaDiff > maxAreaDiff)
                {
                        maxAreaDiff = areaDiff;
                        backObjectsStart = oit + 1;
                }
        }

        // belongs to back bv
        for (oit = objects->begin(); oit != backObjectsStart; ++ oit)
        {
                frontObjects.push_back(*oit);
        }

        // belongs to front bv
        for (oit = backObjectsStart; oit != oit_end; ++ oit)
        {
                backObjects.push_back(*oit);
        }
        
        cout << "front: " << (int)frontObjects.size() << " back: " << (int)backObjects.size() << " " 
                 << backObjects.front()->GetBox().SurfaceArea() - frontObjects.back()->GetBox().SurfaceArea() << endl;
}


inline static float AreaRatio(Intersectable *smallObj, Intersectable *largeObj)
{
        const float areaSmall = smallObj->GetBox().SurfaceArea();
        const float areaLarge = largeObj->GetBox().SurfaceArea();

        return areaSmall / (areaLarge - areaSmall + Limits::Small);
}


bool BvHierarchy::InitialTerminationCriteriaMet(const BvhTraversalData &tData) const
{
        const bool terminationCriteriaMet =
                        (0
                    || ((int)tData.mNode->mObjects.size() < mInitialMinObjects)
                        || (tData.mNode->mObjects.back()->GetBox().SurfaceArea() < mInitialMinArea)
                        || (AreaRatio(tData.mNode->mObjects.front(), tData.mNode->mObjects.back()) > mInitialMaxAreaRatio)
                        );

        cout << "criteria met: "<< terminationCriteriaMet << "\n" 
                 << "size: " << (int)tData.mNode->mObjects.size() << " max: " << mInitialMinObjects << endl
                 << "ratio: " << AreaRatio(tData.mNode->mObjects.front(), tData.mNode->mObjects.back()) << " max: " << mInitialMaxAreaRatio << endl
                 << "area: " << tData.mNode->mObjects.back()->GetBox().SurfaceArea() << " max: " << mInitialMinArea << endl << endl;

        return terminationCriteriaMet;
}


// HACK
float BvHierarchy::GetTriangleSizeIncrementially(BvhNode *node) const
{
        if (node->mRenderCost < 0)
        {
                //cout <<"p";
                if (node->IsLeaf())
                {
                        BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
                        node->mRenderCost = (float)leaf->mObjects.size();
                }
                else
                {
                        BvhInterior *interior = static_cast<BvhInterior *>(node);
                
                        node->mRenderCost = GetTriangleSizeIncrementially(interior->GetFront()) + 
                                                                GetTriangleSizeIncrementially(interior->GetBack());
                }
        }

        return node->mRenderCost;
}


void BvHierarchy::Compress()
{
}


void BvHierarchy::SetUniqueNodeIds()
{
        // export bounding boxes
        vector<BvhNode *> nodes;

        // hack: should also expect interior nodes
        CollectNodes(mRoot, nodes);

        vector<BvhNode *>::const_iterator oit, oit_end = nodes.end();

        int id = 0;

        for (oit = nodes.begin(); oit != oit_end; ++ oit, ++ id)
        {
                (*oit)->SetId(id);
        }
}


}