source: GTP/trunk/Lib/Vis/Preprocessing/src/ViewCell.cpp @ 2071

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


#include "ViewCell.h"
#include "Mesh.h"
#include "Intersectable.h"
#include "KdTree.h"
#include "Triangle3.h"
#include "common.h"
#include "Environment.h"
#include "ViewCellsManager.h"
#include "Exporter.h"
#include "BvHierarchy.h"



namespace GtpVisibilityPreprocessor {



template <typename T> class myless
{
public:
        bool operator() (T v1, T v2) const
        {
                return (v1->GetMergeCost() < v2->GetMergeCost());
        }
};


typedef priority_queue<ViewCell *, vector<ViewCell *>, 
                                           myless<vector<ViewCell *>::value_type> > TraversalQueue;

  int ViewCell::sMailId = 10000;//2147483647;
  int ViewCell::sReservedMailboxes = 1;


float MergeCandidate::sRenderCostWeight = 0;


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

        return (float)pvs;
}

/** Counts contribution of the view cell to the pvs.
*/
inline int CountPvsContribution(ViewCell *vc)
{
        int count = 0;

        ObjectPvsIterator pit = vc->GetPvs().GetIterator();

        while (pit.HasMoreEntries())
        {
                ObjectPvsEntry entry = pit.Next();

                if (!entry.mObject->Mailed())
                {
                        entry.mObject->Mail();
                        ++ count;
                }
        }

        return count;
}


/// Fast computation of merged pvs size
static float ComputeMergedPvsCost(const ObjectPvs &pvs1, 
                                                                  const ObjectPvs &pvs2)
{
        // add first pvs
        float pvs = (float)pvs1.GetSize();

        Intersectable::NewMail();

        // mail all objects in first pvs
        ObjectPvsIterator pit = pvs1.GetIterator();

        while (pit.HasMoreEntries())
        {
                ObjectPvsEntry entry = pit.Next();
                entry.mObject->Mail();
        }

        pit = pvs2.GetIterator();

        while (pit.HasMoreEntries())
        {
                ObjectPvsEntry entry = pit.Next();

                Intersectable *obj = entry.mObject;
                if (!obj->Mailed())
                        ++ pvs;
        }

        return pvs;
}


ViewCell::ViewCell(): 
MeshInstance(NULL), 
mArea(-1),
mVolume(-1),
mValid(true),
mParent(NULL),
mMergeCost(0),
mPvsCost(0),
mEntriesInPvs(0),
mPvsSizeValid(false),
mFilteredPvsSize(0)
{
        mId = -100;
}

ViewCell::ViewCell(Mesh *mesh): 
MeshInstance(mesh), 
mArea(-1),
mVolume(-1),
mValid(true),
mParent(NULL),
mMergeCost(0),
mPvsCost(0),
mPvsSizeValid(false),
mFilteredPvsSize(0)
//mMailbox(0)
{
        mId = -100;
}


ViewCell::~ViewCell()
{
}


const ObjectPvs &ViewCell::GetPvs() const
{
  return mPvs;
}


ObjectPvs &ViewCell::GetPvs()
{
  return mPvs;
}


void ViewCell::SetPvs(const ObjectPvs &pvs)
{
        mPvs = pvs;
}


int ViewCell::Type() const 
{ 
        return VIEW_CELL; 
}


float ViewCell::GetVolume() const
{
        return mVolume;
}


void ViewCell::SetVolume(float volume) 
{
        mVolume = volume;
}


void ViewCell::SetMesh(Mesh *mesh)
{
        mMesh = mesh;
}


float ViewCell::GetArea() const
{
        return mArea;
}


void ViewCell::SetArea(float area) 
{
        mArea = area;
}


void ViewCell::SetColor(const RgbColor &color) 
{
        mColor = color;
}


RgbColor ViewCell::GetColor() const 
{
        return mColor;
}


void ViewCell::SetValid(const bool valid)
{
        mValid = valid;
}


bool ViewCell::GetValid() const
{
        return mValid;
}


void ViewCell::SetParent(ViewCellInterior *parent)
{
        mParent = parent;
}


bool ViewCell::IsRoot() const
{
        return !mParent;
}


ViewCellInterior *ViewCell::GetParent() const
{
        return mParent;
}


void ViewCell::SetMergeCost(const float mergeCost)
{
        mMergeCost = mergeCost;
}


float ViewCell::GetRenderCost() const
{
        return mPvsCost * GetVolume();
}


float ViewCell::GetMergeCost() const
{
        return mMergeCost;
}


bool ViewCell::AddPvsSample(Intersectable *sample,
                                                        const float pdf, 
                                                        float &contribution)
{ 
        const bool result = mPvs.AddSample(sample, pdf);//, contribution);
        // have to recompute pvs size
        mPvsSizeValid = false;

        return result;
}



/************************************************************************/
/*                class ViewCellInterior implementation                 */
/************************************************************************/


ViewCellInterior::ViewCellInterior()
{
}


ViewCellInterior::~ViewCellInterior()
{
        ViewCellContainer::const_iterator it, it_end = mChildren.end();

        for (it = mChildren.begin(); it != it_end; ++ it)
        {
                delete (*it);
        }
}


ViewCellInterior::ViewCellInterior(Mesh *mesh): 
ViewCell(mesh)
{
}


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


void ViewCellInterior::SetupChildLink(ViewCell *vc)
{
    mChildren.push_back(vc);
    vc->SetParent(this);
}


void ViewCellInterior::RemoveChildLink(ViewCell *l)
{
        // erase leaf from old view cell
        ViewCellContainer::iterator it = mChildren.begin();

        for (; (*it) != l; ++ it);
        if (it == mChildren.end())
                Debug << "error" << endl;
        else
                mChildren.erase(it);
}


void ViewCellInterior::ReplaceChildLink(ViewCell *prev, ViewCell *cur)
{
        // erase leaf from old view cell
        ViewCellContainer::iterator it = mChildren.begin();

        for (; (*it) != prev; ++ it);
        if (it == mChildren.end())
        {
                Debug << "error: child link not found" << endl;
        }
        else
        {
                (*it) = cur;
        }
}



/************************************************************************/
/*                class ViewCellsStatistics implementation              */
/************************************************************************/


void ViewCellsStatistics::Print(ostream &app) const
{
        app << "=========== View Cells Statistics ===============\n";

        app << setprecision(4);

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

        app << "#N_OVERALLPVS ( cost of the PVS )\n" << pvsCost << endl;

        //app << "#N_PVSENTRIES ( entries in the PVS)\n" << pvsEntries << endl;

        app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;

        app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;

        app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;

        app << "#N_PEMPTYPVS ( view cells with empty PVS )\n" << emptyPvs << endl;

        app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;

        app << "#N_AVGLEAVES (average number of leaves per view cell )\n" << AvgLeaves() << endl;

        app << "#N_MAXLEAVES ( maximal number of leaves per view cell )\n" << maxLeaves << endl;
        
        app << "#N_INVALID ( number of invalid view cells )\n" << invalid << endl;

        app << "========== End of View Cells Statistics ==========\n";
}


/*************************************************************************/
/*                    class ViewCellsTree implementation                 */
/*************************************************************************/


ViewCellsTree::ViewCellsTree(ViewCellsManager *vcm):
mRoot(NULL),
mUseAreaForPvs(false),
mViewCellsManager(vcm),
#if 0
mViewCellsStorage(PVS_IN_INTERIORS)
#else
mViewCellsStorage(PVS_IN_LEAVES)
#endif
{
        ReadEnvironment();
        MergeCandidate::sRenderCostWeight = mRenderCostWeight;
}


ViewCellsTree::ViewCellsTree():
mRoot(NULL),
mUseAreaForPvs(false),
mViewCellsManager(NULL),
#if 0
mViewCellsStorage(PVS_IN_INTERIORS)
#else
mViewCellsStorage(PVS_IN_LEAVES)
#endif
{
        ReadEnvironment();
        MergeCandidate::sRenderCostWeight = mRenderCostWeight;
}


void ViewCellsTree::ReadEnvironment()
{
        Environment::GetSingleton()->GetBoolValue("ViewCells.Visualization.exportMergedViewCells", mExportMergedViewCells);
        Environment::GetSingleton()->GetFloatValue("ViewCells.maxStaticMemory", mMaxMemory);

        //-- merge options
        Environment::GetSingleton()->GetFloatValue("ViewCells.PostProcess.renderCostWeight", mRenderCostWeight);
        Environment::GetSingleton()->GetIntValue("ViewCells.PostProcess.minViewCells", mMergeMinViewCells);
        Environment::GetSingleton()->GetFloatValue("ViewCells.PostProcess.maxCostRatio", mMergeMaxCostRatio);
        Environment::GetSingleton()->GetBoolValue("ViewCells.PostProcess.refine", mRefineViewCells);    

        Environment::GetSingleton()->GetIntValue("ViewCells.PostProcess.maxMergesPerPass", mMaxMergesPerPass);
        Environment::GetSingleton()->GetFloatValue("ViewCells.PostProcess.avgCostMaxDeviation", mAvgCostMaxDeviation);

        Debug << "============= view cell tree options ================\n";
        Debug << "minimum view cells: " << mMergeMinViewCells << endl;
        Debug << "max cost ratio: " << mMergeMaxCostRatio << endl;
        Debug << "max memory: " << mMaxMemory << endl;
        Debug << "refining view cells: " << mRefineViewCells << endl;
        Debug << "=========== end view cell tree options ===============\n";
}


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


int ViewCellsTree::GetNumInitialViewCells(ViewCell *vc) const
{
        int vcSize = 0;

        stack<ViewCell *> tstack;

        tstack.push(vc);

        while (!tstack.empty())
        {
                ViewCell *vc = tstack.top();
                tstack.pop();

                if (vc->IsLeaf())
                {
                        ++ vcSize;
                }
                else
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
                        
                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }
                        
                }
        }

        return vcSize;
}


void ViewCellsTree::CollectLeaves(ViewCell *vc, ViewCellContainer &leaves) const
{
        stack<ViewCell *> tstack;

        tstack.push(vc);

        while (!tstack.empty())
        {
                ViewCell *vc = tstack.top();
                tstack.pop();

                if (vc->IsLeaf())
                {
                        leaves.push_back(vc);
                }
                else
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }
                }
        }
}


ViewCellsTree::~ViewCellsTree()
{
        DEL_PTR(mRoot);
}


int ViewCellsTree::ConstructMergeTree(const VssRayContainer &rays, 
                                                                          const ObjectContainer &objects)
{
        mNumActiveViewCells = (int)mViewCellsManager->GetViewCells().size();

        float variance = 0;
        float totalPvs = 0;
        float totalRenderCost = 0;

        //////////////////
        //-- compute statistics values of initial view cells

        mViewCellsManager->EvaluateRenderStatistics(totalRenderCost,
                                                                                                mExpectedCost,
                                                                                                mDeviation,
                                                                                                variance,
                                                                                                totalPvs,
                                                                                                mAvgRenderCost);

        /////////////
        //-- fill merge queue
        vector<MergeCandidate> candidates;

        mViewCellsManager->CollectMergeCandidates(rays, candidates);

        while(!candidates.empty())
        {
                MergeCandidate mc = candidates.back();
                candidates.pop_back();
                EvalMergeCost(mc);
                mMergeQueue.push(mc);
        }

        Debug << "************************* merge ***********************************" << endl;  
        Debug << "deviation: " << mDeviation << endl;
        Debug << "avg render cost: " << mAvgRenderCost << endl;
        Debug << "expected cost: " << mExpectedCost << endl;


        ViewCellsManager::PvsStatistics pvsStats;
        mViewCellsManager->GetPvsStatistics(pvsStats);

        //static float expectedValue = pvsStats.avgPvs;
        
        //-- the current view cells are kept in this container
        //-- we start with the current view cells from the view cell manager. 
        //-- The active view cells will change with subsequent merges
        
        // todo: should rather take initial view cells
    ViewCellContainer &activeViewCells = mViewCellsManager->GetViewCells();
        
        
        ViewCell::NewMail();

        MergeStatistics mergeStats;
        mergeStats.Start();
        
        long startTime = GetTime();

        mergeStats.collectTime = TimeDiff(startTime, GetTime());
        mergeStats.candidates = (int)mMergeQueue.size();
        startTime = GetTime();

        // frequency stats are updated
        const int statsOut = 500;
        
        // passes are needed for statistics, because we don't want to record
        // every merge
        int pass = 0;
        int mergedPerPass = 0;
        float realExpectedCost = mExpectedCost;
        float realAvgRenderCost = mAvgRenderCost;
        int realNumActiveViewCells = mNumActiveViewCells;
        
        // maximal ratio of old expected render cost to expected render
        // when the the render queue has to be reset.
        int numMergedViewCells = 0;
                

        cout << "actual merge starts now ... " << endl;

        //-- use priority queue to merge leaf pairs

        while (!mMergeQueue.empty())
        {
                //-- reset merge queue if the ratio of current expected cost / real expected cost
                //   too small or after a given number of merges
                if ((mergedPerPass > mMaxMergesPerPass) ||
                        (mAvgCostMaxDeviation > mAvgRenderCost / realAvgRenderCost))
                {
                        Debug << "************ reset queue *****************\n"
                                  << "ratios: " << mAvgCostMaxDeviation 
                                  << " real avg render cost " << realAvgRenderCost << " average render cost " << mAvgRenderCost
                                  << " merged per pass : " << mergedPerPass << " of maximal " << mMaxMergesPerPass << endl;

                        Debug << "Values before reset: "  
                                  << " erc: " << mExpectedCost 
                                  << " avgrc: " << mAvgRenderCost
                                  << " dev: " << mDeviation << endl;
        
                        // adjust render cost
                        ++ pass;

                        mergedPerPass = 0;
                        mExpectedCost = realExpectedCost;
                        mAvgRenderCost = realAvgRenderCost;
                        mNumActiveViewCells = realNumActiveViewCells;
                        
                        const int numMergedViewCells = UpdateActiveViewCells(activeViewCells);
                
                        /////////////////
                        //-- reset / refine the view cells
                        //-- priorities are recomputed
                        //-- the candidates are put back into merge queue

                        if (mRefineViewCells)
                                RefineViewCells(rays, objects);
                        else
                                ResetMergeQueue();

                        Debug << "Values after reset: "  
                                  << " erc: " << mExpectedCost 
                                  << " avg: " << mAvgRenderCost 
                                  << " dev: " << mDeviation << endl;

                        if (mExportMergedViewCells)
                        {
                                ExportMergedViewCells(activeViewCells, objects, numMergedViewCells);
                        }
                }


                MergeCandidate mc = mMergeQueue.top();
                mMergeQueue.pop();
        
                // both view cells equal because of previous merges
                // NOTE: do I really still need this? probably cannot happen!!
                if (mc.mLeftViewCell == mc.mRightViewCell)
                        continue;

                if (mc.IsValid())
                {
                        ViewCell::NewMail();

                        //-- update statistical values
                        -- realNumActiveViewCells;
                        ++ mergeStats.merged;
                        ++ mergedPerPass;

                        const float renderCostIncr = mc.GetRenderCost();
                        const float mergeCostIncr = mc.GetMergeCost();

                        totalRenderCost += renderCostIncr;
                        mDeviation += mc.GetDeviationIncr();

                                                
                        //-- merge the view cells of leaf1 and leaf2
                        float pvsDiff;
                        ViewCellInterior *mergedVc = 
                                MergeViewCells(mc.mLeftViewCell, mc.mRightViewCell, pvsDiff);   

                        // total render cost and deviation has changed
                        // real expected cost will be larger than expected cost used for the
                        // cost heuristics, but cannot recompute costs on each increase of the 
                        // expected cost
                        totalPvs += pvsDiff;
                        realExpectedCost = totalRenderCost / (float)realNumActiveViewCells;
                        realAvgRenderCost = (float)totalPvs / (float)realNumActiveViewCells;
        
                        // set merge cost to this node for priority traversal
                        mergedVc->SetMergeCost(totalRenderCost);
                        
                        // check if "siblings (back and front node of the same parent)
                        if (0) ++ mergeStats.siblings;

                        // set the cost for rendering a view cell
                        mergedVc->SetCost(realExpectedCost);

                        if ((mergeStats.merged % statsOut) == 0)
                                cout << "merged " << mergeStats.merged << " view cells" << endl;

                }
                else
                { 
                        // merge candidate not valid, because one of the leaves was already
                        // merged with another one => validate and reinsert into queue
                        if (ValidateMergeCandidate(mc))
                        {
                                EvalMergeCost(mc);
                                mMergeQueue.push(mc);
                        }
                }
                
        }

        // adjust stats and reset queue one final time
        mExpectedCost = realExpectedCost;
        mAvgRenderCost = realAvgRenderCost;
        mNumActiveViewCells = realNumActiveViewCells;

        UpdateActiveViewCells(activeViewCells);

        // refine view cells and reset costs
        if (mRefineViewCells)
                RefineViewCells(rays, objects);
        else
                ResetMergeQueue();


        // create a root node if the merge was not done till root level,
        // else take the single node as new root
        if ((int)activeViewCells.size() > 1)
        {
                Debug << "creating root of view cell hierarchy for " 
                          << (int)activeViewCells.size() << " view cells" << endl;
                
                ViewCellInterior *root = mViewCellsManager->MergeViewCells(activeViewCells);
        
                ViewCellContainer::const_iterator it, it_end = root->mChildren.end();

                for (it = root->mChildren.begin(); it != it_end; ++ it)
                        (*it)->SetParent(root);
        
                root->SetMergeCost(totalRenderCost);
                // $$JB keep this 0 temporarilly
                root->SetCost(0.0f);

                mRoot = root;
        }
        // normal case
        else if (!activeViewCells.empty())
        {
                Debug << "setting root of the merge history" << endl;
                mRoot = activeViewCells[0];
        }
        else
        {
                Debug << "big error, root is NULL" << endl;
        }
        
        //-- empty merge queue just in case
        while (!mMergeQueue.empty())
        {
                mMergeQueue.pop();
        }

        // test if computed volumes are correct
        Debug << "volume of the root view cell: " << mRoot->GetVolume() 
                  << " " << mViewCellsManager->GetViewSpaceBox().GetVolume() << endl;
        
        // TODO: delete because makes no sense here
        mergeStats.expectedRenderCost = realExpectedCost;
        mergeStats.deviation = mDeviation;

        // we want to optimize this heuristics
        mergeStats.heuristics = 
                mDeviation * (1.0f - mRenderCostWeight) + 
                mExpectedCost * mRenderCostWeight;

        mergeStats.mergeTime = TimeDiff(startTime, GetTime());
        mergeStats.Stop();
        Debug << mergeStats << endl << endl;

        // assign colors for the view cells so that at least one is always consistent
        AssignRandomColors();

        //TODO: should return sample contributions?
        return mergeStats.merged;
}


ViewCell *ViewCellsTree::GetRoot() const
{
        return mRoot;
}


void ViewCellsTree::ResetMergeQueue()
{
        cout << "reset merge queue ... ";
        
        vector<MergeCandidate> buf;
        buf.reserve(mMergeQueue.size());
                        
        
        // store merge candidates in intermediate buffer
        while (!mMergeQueue.empty())
        {
                MergeCandidate mc = mMergeQueue.top();
                mMergeQueue.pop();
                
                // recalculate cost
                if (ValidateMergeCandidate(mc))
                {
                        EvalMergeCost(mc);
                        buf.push_back(mc);                              
                }
        }

        vector<MergeCandidate>::const_iterator bit, bit_end = buf.end();

        // reinsert back into queue
        for (bit = buf.begin(); bit != bit_end; ++ bit)
        {      
                mMergeQueue.push(*bit); 
        }

        cout << "finished" << endl;
}


float ViewCellsTree::ComputeMergedPvsCost(const ObjectPvs &pvs1, 
                                                                                  const ObjectPvs &pvs2) const
{
        // computes render cost of merge
        float renderCost = 0;

        // compute new pvs size
        Intersectable::NewMail();

        ObjectPvsIterator pit = pvs1.GetIterator();

        // first mail all objects in first pvs
        while (pit.HasMoreEntries())
        {
                ObjectPvsEntry entry = pit.Next();

                Intersectable *obj = entry.mObject;

                obj->Mail();
                renderCost += mViewCellsManager->EvalRenderCost(obj);
        }

        ObjectPvsIterator pit2 = pvs2.GetIterator();

        while (pit2.HasMoreEntries())
        {
                ObjectPvsEntry entry = pit2.Next();
                Intersectable *obj = entry.mObject;

                // test if object already considered    
                if (!obj->Mailed())
                {
                        renderCost += mViewCellsManager->EvalRenderCost(obj);
                }
        }

        return renderCost;
}


int ViewCellsTree::UpdateActiveViewCells(ViewCellContainer &viewCells)
{
        int numMergedViewCells = 0;

        Debug << "updating active vc: " << (int)viewCells.size() << endl;
        
        // find all already merged view cells and remove them from the
        // container view cells
                
        // sort out all view cells which are not active anymore, i.e., they
        // were already part of a merge
        int i = 0;

        ViewCell::NewMail();

        while (1)
        {
                // remove all merged view cells from end of the vector
                while (!viewCells.empty() && (viewCells.back()->GetParent()))
                {
                        viewCells.pop_back();
                }

                // all merged view cells have been found
                if (i >= (int)viewCells.size()) 
                        break;

                // already merged this view cell, put it to end of vector
                if (viewCells[i]->GetParent())
                        swap(viewCells[i], viewCells.back());
                
                // mail view cell that it has not been merged
                viewCells[i]->Mail();

                // increase loop counter
                ++ i;
        }


        // add new view cells to container only if they don't have been 
        // merged in the mean time
        ViewCellContainer::const_iterator ait, ait_end = mMergedViewCells.end();

        for (ait = mMergedViewCells.begin(); ait != ait_end; ++ ait)
        {
                ViewCell *vc = mMergedViewCells.back();
                if (!vc->GetParent() && !vc->Mailed())
                {
                        vc->Mail();
                        viewCells.push_back(vc);
                        ++ numMergedViewCells;
                }
        }

        // dispose old merged view cells
        mMergedViewCells.clear();

        // update standard deviation
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
        
        mDeviation = 0;

        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
        {
                const float lower = (float)mViewCellsManager->GetMinPvsSize();
                const float upper = (float)mViewCellsManager->GetMaxPvsSize();

                const float penalty = EvalPvsPenalty((*vit)->GetPvs().EvalPvsCost(), lower, upper);
                
                mDeviation += fabs(mAvgRenderCost - penalty);
        }

        mDeviation /= (float)viewCells.size();
        
        return numMergedViewCells;
}


void ViewCellsTree::ExportMergedViewCells(ViewCellContainer &viewCells, 
                                                                                  const ObjectContainer &objects,
                                                                                  const int numMergedViewCells)
{
        

        char s[64];

        sprintf(s, "merged_viewcells%07d.x3d", (int)viewCells.size());
        Exporter *exporter = Exporter::GetExporter(s);

        if (exporter)
        {
                cout << "exporting " << (int)viewCells.size() << " merged view cells ... ";
                exporter->ExportGeometry(objects);
                //Debug << "vc size " << (int)viewCells.size() << " merge queue size: " << (int)mMergeQueue.size() << endl;
                ViewCellContainer::const_iterator it, it_end = viewCells.end();

                int i = 0;
                for (it = viewCells.begin(); it != it_end; ++ it)
                {
                        Material m;
                        // assign special material to new view cells
                        // new view cells are on the back of container
                        if (i ++ >= (viewCells.size() - numMergedViewCells))
                        {
                                //m = RandomMaterial();
                                m.mDiffuseColor.r = RandomValue(0.5f, 1.0f);
                                m.mDiffuseColor.g = RandomValue(0.5f, 1.0f);
                                m.mDiffuseColor.b = RandomValue(0.5f, 1.0f);
                        }
                        else
                        {
                                float col = RandomValue(0.1f, 0.4f);
                                m.mDiffuseColor.r = col;
                                m.mDiffuseColor.g = col;
                                m.mDiffuseColor.b = col;
                        }

                        exporter->SetForcedMaterial(m);
                        mViewCellsManager->ExportViewCellGeometry(exporter, *it, NULL, NULL);
                }

                delete exporter;
                cout << "finished" << endl;
        }
}


// TODO: should be done in view cells manager
ViewCellInterior *ViewCellsTree::MergeViewCells(ViewCell *left, 
                                                                                                ViewCell *right, 
                                                                                                float &pvsDiff) //const
{
        // create merged view cell
        ViewCellInterior *vc = 
                mViewCellsManager->MergeViewCells(left, right);

        // if merge was unsuccessful
        if (!vc) return NULL;

        // set to the new parent view cell
        left->SetParent(vc);
        right->SetParent(vc);

        
        if (mUseAreaForPvs)
        {
                // set new area of view cell
                // not not correct, but costly to compute real area!!
                vc->SetArea(left->GetArea() + right->GetArea());
        }
        else
        {       // set new volume of view cell
                vc->SetVolume(left->GetVolume() + right->GetVolume());
        }

        
        // important so other merge candidates sharing this view cell
        // are notified that the merge cost must be updated!!
        vc->Mail();

        const float pvs1 = left->GetPvs().EvalPvsCost();
        const float pvs2 = right->GetPvs().EvalPvsCost();


        // the new view cells are stored in this container
        mMergedViewCells.push_back(vc);

        pvsDiff = vc->GetPvs().EvalPvsCost() - pvs1 - pvs2;


        // don't store pvs in interior cells, just a scalar
        if (mViewCellsStorage == PVS_IN_LEAVES)
        {
                // set scalars
                mViewCellsManager->UpdateScalarPvsSize(left, 
                                left->GetPvs().EvalPvsCost(), 
                                left->GetPvs().GetSize());
                        
                // remove pvs, we don't store interior pvss
                if (!left->IsLeaf())
                {
                        left->GetPvs().Clear();
                }

                // set scalars
                mViewCellsManager->UpdateScalarPvsSize(right, 
                        right->GetPvs().EvalPvsCost(), 
                        right->GetPvs().GetSize());

                right->mPvsSizeValid = true;
                
                // remove pvs, we don't store interior pvss
                if (!right->IsLeaf())
                {
                        right->GetPvs().Clear();
                }
        }

        return vc;
}


int ViewCellsTree::RefineViewCells(const VssRayContainer &rays, 
                                                                   const ObjectContainer &objects)
{
        Debug << "refining " << (int)mMergeQueue.size() << " candidates " << endl;

        // intermediate buffer for shuffled view cells
        vector<MergeCandidate> buf;
        buf.reserve(mMergeQueue.size());
                        
        // Use priority queue of remaining leaf pairs 
        // The candidates either share the same view cells or
        // are border leaves which share a boundary.
        // We test if they can be shuffled, i.e.,
        // either one leaf is made part of one view cell or the other
        // leaf is made part of the other view cell. It is tested if the
        // remaining view cells are "better" than the old ones.
        
        const int numPasses = 3;
        int pass = 0;
        int passShuffled = 0;
        int shuffled = 0;
        int shuffledViewCells = 0;

        ViewCell::NewMail();
        
        while (!mMergeQueue.empty())
        {
                MergeCandidate mc = mMergeQueue.top();
                mMergeQueue.pop();

                // both view cells equal or already shuffled
                if ((mc.GetLeftViewCell() == mc.GetRightViewCell()) ||
                        mc.GetLeftViewCell()->IsLeaf() || mc.GetRightViewCell()->IsLeaf())
                {                       
                        continue;
                }

                // candidate for shuffling
                const bool wasShuffled = ShuffleLeaves(mc);
                
                // shuffled or put into other queue for further refine
                if (wasShuffled)
                {
                        ++ passShuffled;

                        if (!mc.GetLeftViewCell()->Mailed())
                        {
                                mc.GetLeftViewCell()->Mail();
                                ++ shuffledViewCells;
                        }
                        if (!mc.GetRightViewCell()->Mailed())
                        {
                                mc.GetRightViewCell()->Mail();
                                ++ shuffledViewCells;
                        }
                }

                // put back into intermediate vector
                buf.push_back(mc);
        }


        //-- in the end, the candidates must be in the mergequeue again
        //   with the correct cost

        cout << "reset merge queue ... ";
        
        
        vector<MergeCandidate>::const_iterator bit, bit_end = buf.end();
        
        for (bit = buf.begin(); bit != bit_end; ++ bit)
        {    
                MergeCandidate mc = *bit;
                // recalculate cost
                if (ValidateMergeCandidate(mc))
                {
                        EvalMergeCost(mc);
                        mMergeQueue.push(mc);   
                }
        }

        cout << "finished" << endl;

        return shuffledViewCells;
}


inline int AddedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
{
        ObjectPvs interPvs;
        ObjectPvs::Merge(interPvs, pvs1, pvs2);

        return (int)interPvs.GetSize();
}


// recomputes pvs size minus pvs of leaf l
#if 0
inline int SubtractedPvsSize(BspViewCell *vc, BspLeaf *l, const ObjectPvs &pvs2)
{
        ObjectPvs pvs;
        vector<BspLeaf *>::const_iterator it, it_end = vc->mLeaves.end();
        for (it = vc->mLeaves.begin(); it != vc->mLeaves.end(); ++ it)
                if (*it != l)
                        pvs.AddPvs(*(*it)->mPvs);
        return pvs.GetSize();
}
#endif


// computes pvs1 minus pvs2
inline int SubtractedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
{
        return pvs1.SubtractPvs(pvs2);
}


float ViewCellsTree::EvalShuffleCost(ViewCell *leaf, 
                                                                         ViewCellInterior *vc1, 
                                                                         ViewCellInterior *vc2) const
{
        //const int pvs1 = SubtractedPvsSize(vc1, leaf, *leaf->mPvs);
        const float pvs1 = (float)SubtractedPvsSize(vc1->GetPvs(), leaf->GetPvs());
        const float pvs2 = (float)AddedPvsSize(vc2->GetPvs(), leaf->GetPvs());

        const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
        const float upperPvsLimit = (float)mViewCellsManager->GetMaxPvsSize();

        const float pvsPenalty1 = 
                EvalPvsPenalty(pvs1, lowerPvsLimit, upperPvsLimit);

        const float pvsPenalty2 =
                EvalPvsPenalty(pvs2, lowerPvsLimit, upperPvsLimit);


        // don't shuffle leaves with pvs > max
        if (0 && (pvs1 + pvs2 > mViewCellsManager->GetMaxPvsSize()))
        {
                return 1e20f;
        }

        float p1, p2;

    if (mUseAreaForPvs)
        {
                p1 = vc1->GetArea() - leaf->GetArea();
                p2 = vc2->GetArea() + leaf->GetArea();
        }
        else
        {
                p1 = vc1->GetVolume() - leaf->GetVolume();
                p2 = vc2->GetVolume() + leaf->GetVolume();
        }

        const float renderCost1 = pvsPenalty1 * p1;
        const float renderCost2 = pvsPenalty2 * p2;

        float dev1, dev2;

        if (1)
        {
                dev1 = fabs(mAvgRenderCost - pvsPenalty1);
                dev2 = fabs(mAvgRenderCost - pvsPenalty2);
        }
        else
        {
                dev1 = fabs(mExpectedCost - renderCost1);
                dev2 = fabs(mExpectedCost - renderCost2);
        }
        
        return mRenderCostWeight * (renderCost1 + renderCost2) +
                  (1.0f - mRenderCostWeight) * (dev1 + dev2) / (float)mNumActiveViewCells;
}


void ViewCellsTree::ShuffleLeaf(ViewCell *leaf,
                                                                ViewCellInterior *vc1, 
                                                                ViewCellInterior *vc2) const
{
        // compute new pvs and area
        // TODO change
        vc1->GetPvs().SubtractPvs(leaf->GetPvs());
        
        ObjectPvs interPvs;
        ObjectPvs::Merge(interPvs, vc2->GetPvs(), leaf->GetPvs());
        vc2->SetPvs(interPvs);

        if (mUseAreaForPvs)
        {
                vc1->SetArea(vc1->GetArea() - leaf->GetArea());
                vc2->SetArea(vc2->GetArea() + leaf->GetArea());
        }
        else
        {
                vc1->SetVolume(vc1->GetVolume() - leaf->GetVolume());
                vc2->SetVolume(vc2->GetVolume() + leaf->GetVolume());
        }

        
        ViewCellInterior *p = static_cast<ViewCellInterior *>(leaf->GetParent());

        p->RemoveChildLink(leaf);
        vc2->SetupChildLink(leaf);
}


bool ViewCellsTree::ShuffleLeaves(MergeCandidate &mc) const
{
        float cost1, cost2;

        ViewCellInterior *vc1 = static_cast<ViewCellInterior *>(mc.GetLeftViewCell());
        ViewCellInterior *vc2 = static_cast<ViewCellInterior *>(mc.GetRightViewCell());

        ViewCell *leaf1 = mc.GetInitialLeftViewCell();
        ViewCell *leaf2 = mc.GetInitialRightViewCell();

        //-- first test if shuffling would decrease cost
        cost1 = GetCostHeuristics(vc1);
        cost2 = GetCostHeuristics(vc2);

        const float oldCost = cost1 + cost2;
        
        float shuffledCost1 = Limits::Infinity;
        float shuffledCost2 = Limits::Infinity;

        if (leaf1)
                shuffledCost1 = EvalShuffleCost(leaf1, vc1, vc2);       
        if (leaf2)
                shuffledCost2 = EvalShuffleCost(leaf2, vc2, vc1);

        // if cost of shuffle is less than old cost => shuffle
        if ((oldCost <= shuffledCost1) && (oldCost <= shuffledCost2))
                return false;
        
        if (shuffledCost1 < shuffledCost2)
        {
                if (leaf1)
                        ShuffleLeaf(leaf1, vc1, vc2);
                mc.mInitialLeftViewCell = NULL;
        }
        else
        {
                if (leaf2)
                        ShuffleLeaf(leaf2, vc2, vc1);
                mc.mInitialRightViewCell = NULL;
        }

        return true;
}


float ViewCellsTree::GetVariance(ViewCell *vc) const
{
        const float upper = (float)mViewCellsManager->GetMaxPvsSize();
        const float lower = (float)mViewCellsManager->GetMinPvsSize();

        if (1)
        {
                const float penalty = 
                        EvalPvsPenalty(GetPvsCost(vc), lower, upper);

                return (mAvgRenderCost - penalty) * (mAvgRenderCost - penalty) / 
                        (float)mNumActiveViewCells;
        }

    const float leafCost = GetRenderCost(vc);
        return (mExpectedCost - leafCost) * (mExpectedCost - leafCost);
}


float ViewCellsTree::GetDeviation(ViewCell *vc) const
{
        const float upper = (float)mViewCellsManager->GetMaxPvsSize();
        const float lower = (float)mViewCellsManager->GetMinPvsSize();

        if (1)
        {
                const float penalty = EvalPvsPenalty(GetPvsCost(vc), lower, upper);
                return fabs(mAvgRenderCost - penalty) / (float)mNumActiveViewCells;
        }

    const float renderCost = GetRenderCost(vc);
        return fabs(mExpectedCost - renderCost);
}


float ViewCellsTree::GetRenderCost(ViewCell *vc) const
{
        if (mUseAreaForPvs)
        {
                return vc->GetPvs().EvalPvsCost() * vc->GetArea();
        }

        return vc->GetPvs().EvalPvsCost() * vc->GetVolume();
}


float ViewCellsTree::GetCostHeuristics(ViewCell *vc) const
{
        return GetRenderCost(vc) * mRenderCostWeight + 
                   GetDeviation(vc) * (1.0f - mRenderCostWeight);
}


bool ViewCellsTree::ValidateMergeCandidate(MergeCandidate &mc) const
{
        // propagate up so we have only the active view cells
        while (mc.mLeftViewCell->mParent)
        {
                mc.mLeftViewCell = mc.mLeftViewCell->mParent;
        }

        while (mc.mRightViewCell->mParent)
        {
                mc.mRightViewCell = mc.mRightViewCell->mParent;
        }

        // this view cell was already merged
        //return mc.mLeftViewCell && (mc.mLeftViewCell != mc.mRightViewCell);
        return mc.mLeftViewCell != mc.mRightViewCell;
}


void ViewCellsTree::EvalMergeCost(MergeCandidate &mc) const
{
        ///////////////////
        //-- compute pvs difference

        float newPvs;
        
        // not valid if not using const cost per object!!
        newPvs = ComputeMergedPvsCost(mc.mLeftViewCell->GetPvs(), mc.mRightViewCell->GetPvs());


        const float newPenalty = EvalPvsPenalty(newPvs, 
                                                                                        (float)mViewCellsManager->GetMinPvsSize(),
                                                                                        (float)mViewCellsManager->GetMaxPvsSize());

        ViewCell *vc1 = mc.mLeftViewCell;
        ViewCell *vc2 = mc.mRightViewCell;

        //-- compute ratio of old cost
        //-- (i.e., added size of left and right view cell times pvs size)
        //-- to new rendering cost (i.e, size of merged view cell times pvs size)
        const float oldCost = GetRenderCost(vc1) + GetRenderCost(vc2);

    const float newCost = mUseAreaForPvs ? 
                (float)newPenalty * (vc1->GetArea() + vc2->GetArea()) :
                (float)newPenalty * (vc1->GetVolume() + vc2->GetVolume());


        // strong penalty if pvs size too large
        if (0 && (newPvs > mViewCellsManager->GetMaxPvsSize()))
        {
                mc.mRenderCost = 1e20f;
        }
        else
        {
                mc.mRenderCost = (newCost - oldCost) / 
                        mViewCellsManager->GetViewSpaceBox().GetVolume();
        }       
        
        ///////////////////////////
        //-- merge cost also takes deviation into account

        float newDev, oldDev;

        if (1)
                newDev = fabs(mAvgRenderCost - newPenalty) / (float)mNumActiveViewCells;
        else
                newDev = fabs(mExpectedCost - newCost) / (float)mNumActiveViewCells;
        
        oldDev = GetDeviation(vc1) + GetDeviation(vc2);

        // compute deviation increase
        mc.mDeviationIncr = newDev - oldDev;
        
        //Debug << "render cost: " << mc.mRenderCost * mRenderCostWeight << endl;
        //Debug << "standard deviation: " << mc.mDeviationIncr * mRenderCostWeight << endl;
}


void ViewCellsTree::SetViewCellsStorage(int stype)
{
        if (mViewCellsStorage == stype)
                return;

        // TODO
        switch (stype)
        {
        case COMPRESSED:
                CompressViewCellsPvs(mRoot);
                break;
        default:
                break;
        }

        mViewCellsStorage = stype;
}


void ViewCellsTree::CompressViewCellsPvs(ViewCell *root)
{
        if (!root->IsLeaf())
        {
                ViewCellInterior *interior = static_cast<ViewCellInterior *>(root);

        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
                
                // compress child sets first
                for (it = interior->mChildren.begin(); it != it_end; ++ it)
                {
                        CompressViewCellsPvs(*it);
                }

                /*cout << "***************\nbefore pulling up:\n";
                for (it = interior->mChildren.begin(); it != it_end; ++ it)
                {
                        cout << "vc:\n" << (*it)->GetPvs() << endl;
                }*/

                // compress root node
                PullUpVisibility(interior);

                /*cout << "after pulling up:\n";
                cout << "interior:\n" << interior->GetPvs() << endl;

                for (it = interior->mChildren.begin(); it != it_end; ++ it)
                {
                        cout << "vc:\n" << (*it)->GetPvs() << endl;
                }*/
        }
}


void ViewCellsTree::UpdateStats(ofstream &stats,
                                                                const ViewCellsTreeStats &vcStats)
{
         stats << "#Pass\n" << vcStats.mPass << endl
                   << "#Splits\n" << vcStats.mNumViewCells << endl 
                   << "#RenderCostDecrease\n" << vcStats.mRenderCostDecrease << endl // TODO
                   << "#TotalRenderCost\n" << vcStats.mTotalRenderCost << endl
                   << "#CurrentPvs\n" << vcStats.mCurrentPvsCost << endl
                   << "#ExpectedCost\n" << vcStats.mExpectedCost << endl
                   << "#AvgRenderCost\n" << vcStats.mAvgRenderCost << endl
                   << "#Deviation\n" << vcStats.mDeviation << endl
                   << "#TotalPvs\n" << vcStats.mTotalPvsCost << endl
                   << "#TotalEntriesInPvs\n" << vcStats.mEntriesInPvs << endl
                   << "#Memory\n" << vcStats.mMemoryCost << endl
                   << "#PvsSizeDecrease\n" << vcStats.mPvsSizeDecr << endl
                   << "#Volume\n" << vcStats.mVolume << endl
                   << endl;
}


void ViewCellsTree::ExportStats(const string &mergeStats)
{
        TraversalQueue tqueue;
        tqueue.push(mRoot);

        ViewCellsTreeStats vcStats;
        vcStats.Reset();

        //cout << "exporting stats ... " << endl;
        vcStats.mNumViewCells = 1;
        
        const AxisAlignedBox3 box = mViewCellsManager->GetViewSpaceBox();
        const float vol = box.GetVolume();

        const float rootPvs = GetPvsCost(mRoot);
        const int rootEntries = GetPvsEntries(mRoot);
        
        vcStats.mTotalPvsCost = rootPvs;
        vcStats.mTotalRenderCost = rootPvs;
        vcStats.mAvgRenderCost = rootPvs;
        vcStats.mExpectedCost = rootPvs;
        
        vcStats.mEntriesInPvs = rootEntries;

        ofstream stats;
        stats.open(mergeStats.c_str());

        vcStats.mMemoryCost = (float)vcStats.mEntriesInPvs * ObjectPvs::GetEntrySize();

        /////////////
        //-- first view cell

        UpdateStats(stats, vcStats);
                                
                
        //-- go through tree in the order of render cost decrease
        //-- which is the same order as the view cells were merged
        //-- or the reverse order of subdivision for subdivision-only 
        //-- view cell hierarchies.

        while (!tqueue.empty())
        {
                ViewCell *vc = tqueue.top();
                tqueue.pop();

                if (!vc->IsLeaf()) 
                {       
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

                        const float parentCost = GetPvsCost(interior);
                        const int parentPvsEntries = GetPvsEntries(interior);
            const float parentExpCost = (float)parentCost * interior->GetVolume();

                        float childExpCost = 0;
                        float childCost = 0;
                        int childPvsEntries = 0;

                        -- vcStats.mNumViewCells;

                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                ViewCell *vc = *it;

                                const float pvsCost = GetPvsCost(vc);
                                const int pvsEntries = GetPvsEntries(vc);

                                childExpCost += childCost * vc->GetVolume();
                                childCost += pvsCost;
                                childPvsEntries += pvsEntries;

                                tqueue.push(vc);
                                ++ vcStats.mNumViewCells;
                        }

                        // update stats for this view cell
                        const float costDecr = (parentExpCost - childExpCost) / vol;

                        vcStats.mTotalRenderCost -= costDecr;
                        vcStats.mTotalPvsCost += childCost - parentCost;
                        vcStats.mEntriesInPvs += childPvsEntries - parentPvsEntries;

                        vcStats.mExpectedCost = vcStats.mTotalRenderCost / (float)vcStats.mNumViewCells;
                        vcStats.mAvgRenderCost = vcStats.mTotalPvsCost / (float)vcStats.mNumViewCells;

                        vcStats.mMemoryCost = vcStats.mEntriesInPvs * ObjectPvs::GetEntrySize();

                        UpdateStats(stats, vcStats);
                }
        }

        stats.close();
}


void ViewCellsTree::SetRoot(ViewCell *root)
{
        mRoot = root;
}


void ViewCellsTree::CollectBestViewCellSet(ViewCellContainer &viewCells, 
                                                                                   const int numViewCells)
{
        TraversalQueue tqueue;
        tqueue.push(mRoot);
        
        while (!tqueue.empty())
        {
                ViewCell *vc = tqueue.top();
                tqueue.pop();

                // save the view cells if it is a leaf or if enough view cells have already been traversed
                // because of the priority queue, this will be the optimal set of v
                if (vc->IsLeaf() || ((viewCells.size() + tqueue.size() + 1) >= numViewCells))
                {
                        viewCells.push_back(vc);
                }
                else 
                {       
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tqueue.push(*it);
                        }
                }
        }
}


void ViewCellsTree::PullUpVisibility(ViewCellInterior *interior)
{
        const int mail = (int)interior->mChildren.size();
        Intersectable::NewMail(mail);

        ViewCellContainer::const_iterator cit, cit_end = interior->mChildren.end();

    // mail all objects in the leaf sets
        // we are interested in the objects which are present in all leaves
        // => count how often an object is part of a child set
        for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
        {
                ViewCell *vc = *cit;
                ObjectPvsIterator pit = vc->GetPvs().GetIterator();

                while (pit.HasMoreEntries())
                {               
                        const ObjectPvsEntry &entry = pit.Next();
                        Intersectable *obj = entry.mObject;
                        obj->Mail();
                }
        }

        for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
        {
                ViewCell *vc = *cit;

                ObjectPvsIterator pit = vc->GetPvs().GetIterator();

                while (pit.HasMoreEntries())
                {               
                        const ObjectPvsEntry &entry = pit.Next();
                        Intersectable *obj = entry.mObject;

                        obj->IncMail();
                }
        }

        // reset pvs
    interior->GetPvs().Clear(false);
                
        // only the objects which are present in all leaf pvs 
        // should remain in the parent pvs
        // these are the objects which have been mailed in all children
        for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
        {
                ViewCell *vc = *cit;
                
                ObjectPvsIterator pit = vc->GetPvs().GetIterator();

                while (pit.HasMoreEntries())
                {               
                        ObjectPvsEntry entry = pit.Next();

                        if (entry.mObject->Mailed(mail))
                        {       
                                interior->GetPvs().AddSample(entry.mObject, entry.mData.mSumPdf);
                        }
                }
        }

        // delete entries which are pulled up
        // note: could be inefficent, rather gather unmailed pvs and reassign
        for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
        {
                ViewCell *vc = *cit;
                ObjectPvsIterator pit = vc->GetPvs().GetIterator();

                ObjectPvs newPvs;

                while (pit.HasMoreEntries())
                {               
                        const ObjectPvsEntry &entry = pit.Next();
                        Intersectable *obj = entry.mObject;

                        /*if (obj->Mailed(mail))
                        {
                                vc->GetPvs().RemoveSample(entry.mObject, Limits::Infinity))
                        }*/
                        if (!obj->Mailed(mail))
                        {
                                newPvs.AddSampleDirty(obj, entry.mData.mSumPdf);
                        }
                }

                vc->SetPvs(newPvs);
        }
}


// TODO matt: implement this function for different storing methods
void ViewCellsTree::GetPvs(ViewCell *vc, ObjectPvs &pvs) const
{
        // pvs is stored in each cell => just return pvs
        if (mViewCellsStorage == PVS_IN_INTERIORS)
        {
                pvs = vc->GetPvs();
                return;
        }

        ////////////////
        //-- pvs is not stored with the interiors => reconstruct
        ViewCell *root = vc;
        
        // add pvs from this view cell to root
        while (root->GetParent())
        {
                root = root->GetParent();
                pvs.MergeInPlace(root->GetPvs());
        }

        // add pvs from leaves
        stack<ViewCell *> tstack;
        tstack.push(vc);

        while (!tstack.empty())
        {
                vc = tstack.top();
                tstack.pop();
        
                // add newly found pvs to merged pvs
                pvs.MergeInPlace(vc->GetPvs());
                
                if (!vc->IsLeaf()) // interior cells: go down to leaf level
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }               
                }
        }
}


float ViewCellsTree::GetPvsCostForLeafStorage(ViewCell *vc) const
{
        // pvs is always stored directly in leaves
        if (vc->IsLeaf())
        {
                return vc->GetPvs().EvalPvsCost();
        }
        
        // interior nodes: pvs is either stored as a 
        // scalar or has to be reconstructed from the leaves
        // the stored pvs size is the valid pvs size => just return scalar
        if (vc->mPvsSizeValid)
        {
                return vc->mPvsCost;
        }
        
        // if no valid pvs size stored as a scalar => 
        // compute current pvs size of interior from it's leaf nodes
        ViewCellContainer leaves;
        CollectLeaves(vc, leaves);

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

        ObjectPvs newPvs;

        // sum up uniquely found intersectables
        for (it = leaves.begin(); it != it_end; ++ it)
        {
                newPvs.MergeInPlace((*it)->GetPvs());
        }

        return newPvs.EvalPvsCost();
}


int ViewCellsTree::GetEntriesInPvsForLeafStorage(ViewCell *vc) const
{
        // pvs is always stored directly in leaves
        if (vc->IsLeaf())
        {
                return vc->GetPvs().GetSize();
        }
        
        // interior node

        // interior nodes: pvs is either stored as a scalar or
        // has to be reconstructed from the leaves

        // the stored pvs size is the valid pvs size => just return scalar
        if (vc->mPvsSizeValid)
        {
                return vc->mEntriesInPvs;
        }
        
        int pvsSize = 0;

        // if no valid pvs size stored as a scalar => 
        // compute current pvs size of interior from it's leaf nodes
        ViewCellContainer leaves;
        CollectLeaves(vc, leaves);

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

        // sum different intersectables
        for (it = leaves.begin(); it != it_end; ++ it)
        {
                ObjectPvsIterator oit = (*it)->GetPvs().GetIterator();

                while (oit.HasMoreEntries())
                {
                        ObjectPvsEntry entry = oit.Next();
                        Intersectable *intersect = entry.mObject;

                        if (!intersect->Mailed())
                        {
                                intersect->Mail();
                                ++ pvsSize;
                        }
                }
        }

        return pvsSize;
}


float ViewCellsTree::GetPvsCostForCompressedStorage(ViewCell *vc) const
{
        float pvsCost = 0;

        /////////////
        //-- compressed pvs

        // the stored pvs size is the valid pvs size => just return scalar
        if (vc->mPvsSizeValid)
        {
                pvsCost = vc->mPvsCost;
        }

        // if no pvs size stored, compute new one
        ViewCell *root = vc;
        
        // also add pvs from this view cell to root
        while (root->GetParent())
        {
                root = root->GetParent();
        
                // matt: bug! must evaluate kd pvss also
                pvsCost += CountPvsContribution(root);
        }

        stack<ViewCell *> tstack;
        tstack.push(vc);

        while (!tstack.empty())
        {
                vc = tstack.top();
                tstack.pop();

                // matt: bug! must evaluate kd pvss also
                pvsCost += CountPvsContribution(vc);

                if (!vc->IsLeaf())
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }               
                }
        }

        return pvsCost;
}


int ViewCellsTree::GetEntriesInPvsForCompressedStorage(ViewCell *vc) const
{
        int pvsSize = 0;

        /////////////////
        //-- compressed pvs

        // the stored pvs size is the valid pvs size => just return scalar
        if (vc->mPvsSizeValid)
        {
                pvsSize = vc->mEntriesInPvs;
        }

        // if no pvs size stored, compute new one
        ViewCell *root = vc;
        
        // also add pvs from this view cell to root
        while (root->GetParent())
        {
                root = root->GetParent();
                // count the pvs entries different from the already found ones  
                pvsSize += CountPvsContribution(root);
        }

        stack<ViewCell *> tstack;
        tstack.push(vc);

        while (!tstack.empty())
        {
                vc = tstack.top();
                tstack.pop();
        
                // count the pvs entries different from the already found ones  
                pvsSize += CountPvsContribution(vc);

                if (!vc->IsLeaf())
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }               
                }
        }

        return pvsSize;
}


float ViewCellsTree::GetPvsCost(ViewCell *vc) const
{
        float pvsCost = 0;
        Intersectable::NewMail();

        ////////////////////////////////////////////////
        // for interiors, pvs can be stored using different methods
        //

        switch (mViewCellsStorage)
        {
        case PVS_IN_LEAVES: 
                // pvs is stored only in leaves
                pvsCost = GetPvsCostForLeafStorage(vc);                 
                break;
        case COMPRESSED:
                pvsCost = GetPvsCostForCompressedStorage(vc);
                break;
        case PVS_IN_INTERIORS:
        default:
                // pvs is stored consistently in the tree up
                //to the root just return pvs size 
                pvsCost = vc->GetPvs().EvalPvsCost();   
                break;
        }

        return pvsCost;  
}


int ViewCellsTree::GetPvsEntries(ViewCell *vc) const
{
        int pvsSize = 0;

        Intersectable::NewMail();

        ////////////////////////////////////////////////
        // for interiors, pvs can be stored using different methods
        
        switch (mViewCellsStorage)
        {
        case PVS_IN_LEAVES:
                {
                        //-- store pvs only in leaves
                        pvsSize = GetEntriesInPvsForLeafStorage(vc);
                        break;
                }
        case COMPRESSED:
                {
                        pvsSize = GetEntriesInPvsForCompressedStorage(vc);
                        break;
                }
        case PVS_IN_INTERIORS:
        default:
                // pvs is stored consistently in the tree up to the root
                // just return pvs size 
                pvsSize = vc->GetPvs().GetSize();       
                break;
        }

        return pvsSize;  
}


float ViewCellsTree::GetMemoryCost(ViewCell *vc) const
{
        return (float)CountStoredPvsEntries(vc) * ObjectPvs::GetEntrySize();
}

//#if HAS_TO_BE_REDONE
int ViewCellsTree::CountStoredPvsEntries(ViewCell *root) const
{
        int pvsSize = root->GetPvs().GetSize();

        // recursivly count leaves
        if (!root->IsLeaf())
        {
                ViewCellInterior *interior = static_cast<ViewCellInterior *>(root);

                ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                for (it = interior->mChildren.begin(); it != it_end; ++ it)
                {
                        pvsSize += CountStoredPvsEntries(*it);
                }
        }

        return pvsSize;         
}
//#endif

int ViewCellsTree::ViewCellsStorage() const
{
        return mViewCellsStorage;
}


ViewCell *ViewCellsTree::GetActiveViewCell(ViewCellLeaf *vc) const
{
        return vc->GetActiveViewCell();
}


void ViewCellsTree::PropagatePvs(ViewCell *root)
{       
        ViewCell *viewCell = root;

        // propagate pvs up
        while (viewCell->GetParent())
        {
                ObjectPvs mergedPvs;
                viewCell->GetParent()->GetPvs().MergeInPlace(root->GetPvs());

                viewCell = viewCell->GetParent();
        }

        if (root->IsLeaf())
                return;

        // propagate pvs to the leaves
        stack<ViewCell *> tstack;
        tstack.push(root);

        while (!tstack.empty())
        {
                ViewCell *viewCell = tstack.top();
                tstack.pop();

                if (viewCell != root)
                {
                        viewCell->GetPvs().MergeInPlace(root->GetPvs());
                }

                if (!viewCell->IsLeaf())
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(viewCell);

                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }
                }
        }
}


void ViewCellsTree::AssignRandomColors()
{
        TraversalQueue tqueue;
        tqueue.push(mRoot);
        
        mRoot->SetColor(RandomColor(0.3f, 1.0f));
        
        while (!tqueue.empty())
        {
                ViewCell *vc = tqueue.top();
                tqueue.pop();

                // save the view cells if it is a leaf or if enough view cells 
                // have already been traversed because of the priority queue, 
                // this will be the optimal set of v
                if (!vc->IsLeaf()) 
                {       
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
                  
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
                  
                        float maxProbability = -1.0f;
                  
                        ViewCell *maxViewCell = NULL;
                        for (it = interior->mChildren.begin(); it != it_end; ++ it) 
                        {
                                ViewCell *v = *it;
                          
                                // set random color
                                v->SetColor(RandomColor(0.3f, 1.0f));
                          
                                if (v->GetVolume() > maxProbability) 
                                {
                                        maxProbability = v->GetVolume();
                                        maxViewCell = v;
                                }

                                if (maxViewCell)
                                {
                                        maxViewCell->SetColor(vc->GetColor());
                                }
                                
                                tqueue.push(v);
                        }
                }        
        }
}


/** Get costs resulting from each merge step. 
*/
void ViewCellsTree::GetCostFunction(vector<float> &costFunction)
{
        TraversalQueue tqueue;
        tqueue.push(mRoot);

        int numViewCells = 1;

        const float vol = mViewCellsManager->GetViewSpaceBox().GetVolume();
        const float rootPvs = GetPvsCost(mRoot);
        
        float totalRenderCost;

        float totalPvsCost = rootPvs;
        totalRenderCost = (float)rootPvs;

        costFunction.push_back(totalRenderCost);

        //-- go through tree in the order of render cost decrease
        //-- which is the same order as the view cells were merged
        //-- or the reverse order of subdivision for subdivision-only 
        //-- view cell hierarchies.

        while (!tqueue.empty())
        {
                ViewCell *vc = tqueue.top();
                tqueue.pop();

                if (!vc->IsLeaf()) 
                {       
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

                        const float parentCost = GetPvsCost(interior);
                        const float parentExpCost = parentCost * interior->GetVolume();

                        -- numViewCells;

                        float childExpCost = 0;
                        float childCost = 0;
                        
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                ViewCell *vc = *it;

                                const float pvsCost = GetPvsCost(vc);
                                
                                childExpCost += (float) pvsCost * vc->GetVolume();
                                childCost += pvsCost;
                                
                                tqueue.push(vc);
                                ++ numViewCells;
                        }

                        // update stats for this view cell
                        const float costDecr = (parentExpCost - childExpCost) / vol;
                        totalRenderCost -= costDecr;
                        
                        costFunction.push_back(totalRenderCost);
                }
        }
}


/** Get storage costs resulting from each merge step. 
*/
void ViewCellsTree::GetStorageFunction(vector<int> &storageFunction)
{
        TraversalQueue tqueue;
        tqueue.push(mRoot);

        int numViewCells = 1;

        const float vol = mViewCellsManager->GetViewSpaceBox().GetVolume();
        const int rootEntries = GetPvsEntries(mRoot);

        int entriesInPvs = rootEntries;
    const int entryStorage = sizeof(PvsData) + sizeof(int); // one entry into the pvs

        storageFunction.push_back(rootEntries);

        ////////////
        //-- go through tree in the order of render cost decrease
        //-- which is the same order as the view cells were merged
        //-- or the reverse order of subdivision for subdivision-only 
        //-- view cell hierarchies.

        while (!tqueue.empty())
        {
                ViewCell *vc = tqueue.top();
                tqueue.pop();

                if (!vc->IsLeaf()) 
                {       
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);

                        const float parentPvsCost = GetPvsCost(interior);
                        const int parentPvsEntries = GetPvsEntries(interior);
            const float parentExpCost = (float)parentPvsCost * interior->GetVolume();

                        float childExpCost = 0;
                        int childPvs = 0;
                        int childPvsEntries = 0;

                        -- numViewCells;

                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                ViewCell *vc = *it;

                                const int pvsEntries = GetPvsEntries(vc);
                                childPvsEntries += pvsEntries;

                                tqueue.push(vc);
                                ++ numViewCells;
                        }

                        // update stats for this view cell
                        const float costDecr = (parentExpCost - childExpCost) / vol;

                        entriesInPvs += childPvsEntries - parentPvsEntries;

                        const int storageCost = entriesInPvs * entryStorage;
                        storageFunction.push_back(storageCost);
                }
        }
}



void ViewCellsTree::UpdateViewCellsStats(ViewCell *vc, 
                                                                                 ViewCellsStatistics &vcStat)
{
        ++ vcStat.viewCells;
                
        const float pvsCost = GetPvsCost(vc);

        vcStat.pvsCost += pvsCost;

        if (pvsCost < Limits::Small)
                ++ vcStat.emptyPvs;

        if (pvsCost > vcStat.maxPvs)
                vcStat.maxPvs = pvsCost;

        if (pvsCost < vcStat.minPvs)
                vcStat.minPvs = pvsCost;

        if (!vc->GetValid())
                ++ vcStat.invalid;
}


bool ViewCellsTree::Export(OUT_STREAM &stream, const bool exportPvs)
{
        // export recursivly all view cells from the root
        ExportViewCell(mRoot, stream, exportPvs);

        return true;
}


void ViewCellsTree::CreateUniqueViewCellsIds()
{
        stack<ViewCell *> tstack;

        int currentId = 0;

        tstack.push(mRoot);

        while (!tstack.empty())
        {
                ViewCell *vc = tstack.top();
                tstack.pop();

                if (vc->GetId() != -1) // out of bounds
                        vc->SetId(currentId ++);

                if (!vc->IsLeaf())
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
                        
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();
                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }
                }
        }
}


void ViewCellsTree::ExportPvs(ViewCell *viewCell, OUT_STREAM &stream)
{
        ObjectPvsIterator it = viewCell->GetPvs().GetIterator();

        while (it.HasMoreEntries())
        {
                ObjectPvsEntry entry = it.Next();
                Intersectable *obj = entry.mObject;

                // hack: just output full pvs
                if (0 && (obj->Type() == Intersectable::BVH_INTERSECTABLE))
                {
                        ObjectContainer objects; 
                        BvhNode *node = static_cast<BvhNode *>(obj);
                        node->CollectObjects(objects);
                
                        ObjectContainer::const_iterator oit, oit_end = objects.end();
                        for (oit = objects.begin(); oit != oit_end; ++ oit)
                        {
                                stream << (*oit)->GetId() << " ";
                        }
                }
                else
                {
                        stream << entry.mObject->GetId() << " ";
                }
        }
}


void ViewCellsTree::ExportViewCell(ViewCell *viewCell, 
                                                                   OUT_STREAM &stream, 
                                                                   const bool exportPvs)
{
        if (viewCell->IsLeaf())
        {
                stream << "<Leaf ";
                stream << "id=\"" << viewCell->GetId() << "\" ";
                stream << "active=\"" << static_cast<ViewCellLeaf *>(viewCell)->GetActiveViewCell()->GetId() << "\" ";
                stream << "mergecost=\"" << viewCell->GetMergeCost() << "\" ";
                stream << "pvs=\"";
                
                //-- export pvs, i.e., the ids of the objects in the pvs
                if (exportPvs)
                {
                        ExportPvs(viewCell, stream);
                }
                stream << "\" />" << endl;
        }
        else
        {
                ViewCellInterior *interior = static_cast<ViewCellInterior *>(viewCell);
        
                stream << "<Interior ";
                stream << "id=\"" << viewCell->GetId() << "\" ";
                stream << "mergecost=\"" << viewCell->GetMergeCost() << "\" ";
                stream << "pvs=\"";

                // NOTE: do not export pvss for interior view cells because
                // they can be completely reconstructed from the leaf pvss
                // on the other hand: we could store a tag with the compression scheme,
        // then some scheme were pvs is in the interiors could be used
                if (exportPvs)
                {
                        ExportPvs(viewCell, stream);
                }

                stream << "\" >" << endl;

                //-- recursivly export child view cells
                ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                for (it = interior->mChildren.begin(); it != it_end; ++ it)
                {
                        ExportViewCell(*it, stream, exportPvs);
                }

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


void ViewCellsTree::ResetPvs()
{
        stack<ViewCell *> tstack;

        tstack.push(mRoot);

        while (!tstack.empty())
        {
                ViewCell *vc = tstack.top();
                tstack.pop();

                vc->GetPvs().Clear();
                
                if (!vc->IsLeaf())
                {
                        ViewCellInterior *interior = static_cast<ViewCellInterior *>(vc);
                        ViewCellContainer::const_iterator it, it_end = interior->mChildren.end();

                        for (it = interior->mChildren.begin(); it != it_end; ++ it)
                        {
                                tstack.push(*it);
                        }
                }
        }
}


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


void ViewCellsTree::SetActiveSetToLeaves()
{
        // todo
}



/**************************************************************************/
/*                     MergeCandidate implementation                      */
/**************************************************************************/


MergeCandidate::MergeCandidate(ViewCell *l, ViewCell *r):
mRenderCost(0),
mDeviationIncr(0),
mLeftViewCell(l),
mRightViewCell(r),
mInitialLeftViewCell(l),
mInitialRightViewCell(r)
{
        //EvalMergeCost();
}


void MergeCandidate::SetRightViewCell(ViewCell *v)
{
        mRightViewCell = v;
}


void MergeCandidate::SetLeftViewCell(ViewCell *v)
{
        mLeftViewCell = v;
}


ViewCell *MergeCandidate::GetRightViewCell() const
{
        return mRightViewCell;
}


ViewCell *MergeCandidate::GetLeftViewCell() const
{
        return mLeftViewCell;
}


ViewCell *MergeCandidate::GetInitialRightViewCell() const
{
        return mInitialRightViewCell;
}


ViewCell *MergeCandidate::GetInitialLeftViewCell() const
{
        return mInitialLeftViewCell;
}


bool MergeCandidate::IsValid() const
{
        // if one has a parent, it was already merged
        return !(mLeftViewCell->GetParent() || mRightViewCell->GetParent());
}


float MergeCandidate::GetRenderCost() const
{
        return mRenderCost;
}


float MergeCandidate::GetDeviationIncr() const
{
        return mDeviationIncr;
}


float MergeCandidate::GetMergeCost() const
{
        return mRenderCost * sRenderCostWeight + 
                   mDeviationIncr * (1.0f - sRenderCostWeight);
}



/************************************************************************/
/*                    MergeStatistics implementation                    */
/************************************************************************/


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

        app << setprecision(4);

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

        app << "#N_CCTIME ( Collect candidates time [s] )\n" << collectTime * 1e-3f << " \n";

        app << "#N_MTIME ( Merge time [s] )\n" << mergeTime * 1e-3f << " \n";

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

        app << "#N_CANDIDATES ( Number of merge candidates )\n" << candidates << "\n";

        app << "#N_MERGEDSIBLINGS ( Number of merged siblings )\n" << siblings << "\n";

        app << "#OVERALLCOST ( overall merge cost )\n" << overallCost << "\n";

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

        app << "#MAX_TREEDIST ( Maximal distance in tree of merged leaves )\n" << maxTreeDist << "\n";

        app << "#AVG_TREEDIST ( Average distance in tree of merged leaves )\n" << AvgTreeDist() << "\n";

        app << "#EXPECTEDCOST ( expected render cost )\n" << expectedRenderCost << "\n";

        app << "#DEVIATION ( deviation )\n" << deviation << "\n";

        app << "#HEURISTICS ( heuristics )\n" << heuristics << "\n";
        

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


}