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

#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"


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 = 1;//2147483647;
int ViewCell::sReservedMailboxes = 1;


float MergeCandidate::sRenderCostWeight = 0;


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

        return (float)pvs;
}

/// Counts differences between pvss.
inline int CountDiffPvs(ViewCell *vc)
{
        int count = 0;

        ObjectPvsMap::const_iterator it, it_end = vc->GetPvs().mEntries.end();
        for (it = vc->GetPvs().mEntries.begin(); it != it_end; ++ it)
        {
                if (!(*it).first->Mailed())
                {
                        (*it).first->Mail();
                        ++ count;
                }
        }

        return count;
}


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

        ObjectPvsMap::const_iterator it, it_end =  pvs1.mEntries.end();

        Intersectable::NewMail();

        // mail all objects in first pvs
        for (it = pvs1.mEntries.begin(); it != it_end; ++ it)
        {
                (*it).first->Mail();
        }

        it_end = pvs2.mEntries.end();

        // look if the entries are also in second pvs
        for (it = pvs2.mEntries.begin(); it != it_end; ++ it)
        {
                Intersectable *obj = (*it).first;
                if (!obj->Mailed())
                        ++ pvs;
        }

        return pvs;
}


ViewCell::ViewCell(): 
MeshInstance(NULL), 
mArea(-1),
mVolume(-1),
mValid(true),
mParent(NULL),
mMergeCost(0),
mPvsSize(0),
mEntriesInPvs(0),
mPvsSizeValid(false)
//mMailbox(0)
{
}

ViewCell::ViewCell(Mesh *mesh): 
MeshInstance(mesh), 
mArea(-1),
mVolume(-1),
mValid(true),
mParent(NULL),
mMergeCost(0),
mPvsSize(0),
mPvsSizeValid(false)
//mMailbox(0)
{
}


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 (float)mPvs.GetSize() * GetVolume();
        return (float)mPvsSize * 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);
        mPvsSizeValid = false; // have to recompute pvs size

        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 ( objects in PVS )\n" << pvsSize << 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 = dynamic_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 = dynamic_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;
        int 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);
                
                        
                        //-- resets / refines 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
                        int 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);
                        
                        // HACK
                        //mergedVc->SetMergeCost(1.0f / (float)realNumActiveViewCells);

                        // check if "siblings (back and front node of the same parent)
                        if (0 && mViewCellsManager->EqualToSpatialNode(mergedVc))
                                ++ mergeStats.siblings;
                        // set the coŽst 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 voluje is equal
        Debug << "volume of the root view cell: " << mRoot->GetVolume() << " " << mViewCellsManager->GetViewSpaceBox().GetVolume() << endl;

        //hack!!
        //mRoot->GetPvs().Clear();
        
        // 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
        ObjectPvsMap::const_iterator it, it_end =  pvs1.mEntries.end();

        Intersectable::NewMail();

        // first mail all objects in first pvs
        for (it = pvs1.mEntries.begin(); it != it_end; ++ it)
        {
                Intersectable *obj = (*it).first;

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

        it_end = pvs2.mEntries.end();

        
        for (it = pvs2.mEntries.begin(); it != it_end; ++ it)
        {
                Intersectable *obj = (*it).first;

                // 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 int lower = mViewCellsManager->GetMinPvsSize();
                const int upper = mViewCellsManager->GetMaxPvsSize();

                const float penalty = EvalPvsPenalty((*vit)->GetPvs().CountObjectsInPvs(), 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, 
                                                                                                int &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 int pvs1 = left->GetPvs().CountObjectsInPvs();
        const int pvs2 = right->GetPvs().CountObjectsInPvs();


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

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


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

                // set scalars
                mViewCellsManager->UpdateScalarPvsSize(right, 
                        right->GetPvs().CountObjectsInPvs(), 
                        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)
{
        return pvs1.AddPvs(pvs2);
}


// 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 int pvs1 = SubtractedPvsSize(vc1->GetPvs(), leaf->GetPvs());
        const int pvs2 = AddedPvsSize(vc2->GetPvs(), leaf->GetPvs());

        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
        const int upperPvsLimit = 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());
        vc2->GetPvs().AddPvs(leaf->GetPvs());
        
        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 = dynamic_cast<ViewCellInterior *>(leaf->GetParent());

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


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

        ViewCellInterior *vc1 = dynamic_cast<ViewCellInterior *>(mc.GetLeftViewCell());
        ViewCellInterior *vc2 = dynamic_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 int upper = mViewCellsManager->GetMaxPvsSize();
        const int lower = mViewCellsManager->GetMinPvsSize();

        if (1)
        {
                const float penalty = 
                        EvalPvsPenalty(vc->GetPvs().GetSize(), 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 int upper = mViewCellsManager->GetMaxPvsSize();
        const int lower = mViewCellsManager->GetMinPvsSize();

        if (1)
        {
                const float penalty = EvalPvsPenalty(vc->GetPvs().CountObjectsInPvs(), 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().CountObjectsInPvs() * vc->GetArea();
        }

        return vc->GetPvs().CountObjectsInPvs() * 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
        int newPvs;
        if (1) // not valid if not using const cost per object!!
                newPvs = ComputeMergedPvsSize(mc.mLeftViewCell->GetPvs(), mc.mRightViewCell->GetPvs());
        else
                newPvs = (int)ComputeMergedPvsCost(mc.mLeftViewCell->GetPvs(), mc.mRightViewCell->GetPvs());

        const float newPenalty = EvalPvsPenalty(newPvs, 
                                                                                        mViewCellsManager->GetMinPvsSize(),
                                                                                        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 = dynamic_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);
                }

                // compress root node
                PullUpVisibility(interior);
        }
}


void ViewCellsTree::UpdateStats(ofstream &stats,
                                                                const int pass,
                                                                const int viewCells,
                                                                const float renderCostDecrease,
                                                                const float totalRenderCost,
                                                                const int currentPvs,
                                                                const float expectedCost,
                                                                const float avgRenderCost,
                                                                const float deviation,
                                                                const int totalPvs,
                                                                const int entriesInPvs,
                                                                const int pvsSizeDecr,
                                                                const float volume)
{
         stats << "#Pass\n" << pass << endl
                << "#ViewCells\n" << viewCells << endl 
        << "#RenderCostDecrease\n" << renderCostDecrease << endl // TODO
                << "#TotalRenderCost\n" << totalRenderCost << endl
                << "#CurrentPvs\n" << currentPvs << endl
                << "#ExpectedCost\n" << expectedCost << endl
                << "#AvgRenderCost\n" << avgRenderCost << endl
                << "#Deviation\n" << deviation << endl
                << "#TotalPvs\n" << totalPvs << endl
                << "#TotalEntriesInPvs\n" << entriesInPvs << endl
                << "#PvsSizeDecrease\n" << pvsSizeDecr << endl // TODO
                << "#Volume\n" << volume << endl
                << endl;
}


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

        tqueue.push(mRoot);
        int numViewCells = 1;
        
        const AxisAlignedBox3 box = mViewCellsManager->GetViewSpaceBox();
        const float vol = box.GetVolume();

        const int rootPvs = GetPvsSize(mRoot);
        const int rootEntries = GetPvsEntries(mRoot);

        Debug << "******** Export stats **********" << endl;
        /*Debug << "vsb volume: " << vol << endl;
        Debug << "root volume: " << mRoot->GetVolume() << endl;
        Debug << "root pvs: " << rootPvs << endl;
        */

        float totalRenderCost, avgRenderCost, expectedCost;

        float deviation = 0;
        int totalPvs = rootPvs;
        int entriesInPvs = rootEntries;

        totalRenderCost = avgRenderCost = expectedCost = (float)rootPvs;

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

        //-- first view cell
        UpdateStats(stats,
                                0, numViewCells, 0, totalRenderCost, 
                                rootPvs, expectedCost, avgRenderCost, deviation,
                                totalPvs, entriesInPvs, 0, mRoot->GetVolume());
                

        //-- 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 = dynamic_cast<ViewCellInterior *>(vc);

                        const int parentPvs = GetPvsSize(interior);
                        const int parentPvsEntries = GetPvsEntries(interior);
            const float parentCost = (float)parentPvs * interior->GetVolume();

                        float childCost = 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 pvsSize = GetPvsSize(vc);
                                const int pvsEntries = GetPvsEntries(vc);

                                childCost += (float) pvsSize * vc->GetVolume();
                                childPvs += pvsSize;
                                childPvsEntries += pvsEntries;

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


                        const float costDecr = (parentCost - childCost) / vol;

                        totalRenderCost -= costDecr;
                        totalPvs += childPvs - parentPvs;
                        entriesInPvs += childPvsEntries - parentPvsEntries;

                        expectedCost = totalRenderCost / (float)numViewCells;
                        avgRenderCost = (float)totalPvs / (float)numViewCells;

                        UpdateStats(stats,
                                                0, numViewCells, costDecr, totalRenderCost,
                                                parentPvs, expectedCost, avgRenderCost, deviation,
                        totalPvs, entriesInPvs, childPvs - parentPvs,
                                                vc->GetVolume());

                }
        }

        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 = dynamic_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)
{
        Intersectable::NewMail((int)interior->mChildren.size());

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

        ObjectPvsMap::const_iterator oit;

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

                ObjectPvsMap::const_iterator oit_end = vc->GetPvs().mEntries.end();

                for (oit = vc->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
                {
                        Intersectable *obj = (*oit).first;
                        if ((cit == interior->mChildren.begin()) && !obj->Mailed())
                                obj->Mail();
                        
                        int incm = obj->IncMail();
                }
        }

        interior->GetPvs().Clear();
        
        
        // 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;

                ObjectPvsMap::const_iterator oit_end = vc->GetPvs().mEntries.end();

                for (oit = vc->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
                {               
                        if ((*oit).first->Mailed((int)interior->mChildren.size()))
                        {       
                                interior->GetPvs().AddSample((*oit).first, (*oit).second.mSumPdf);
                        }
                }
        }


        // delete all the objects from the leaf sets which were moved to parent pvs
        ObjectPvsMap::const_iterator oit_end = interior->GetPvs().mEntries.end();

        for (oit = interior->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
        {
                for (cit = interior->mChildren.begin(); cit != cit_end; ++ cit)
                {
                        if (!(*cit)->GetPvs().RemoveSample((*oit).first, Limits::Infinity))
                        {
                                Debug << "should not come here!" << endl;
                        }
                }
        }
}


// 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
        Intersectable::NewMail();

        int pvsSize = 0;
        ViewCell *root = vc;
        
        // add pvs from this view cell to root
        while (root->GetParent())
        {
                root = root->GetParent();
                pvs.AddPvs(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.AddPvs(vc->GetPvs());

                if (!vc->IsLeaf()) // interior cells: go down to leaf level
                {
                        ViewCellInterior *interior = dynamic_cast<ViewCellInterior *>(vc);

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

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


int ViewCellsTree::GetPvsSizeForLeafStorage(ViewCell *vc) const
{
        // pvs is always stored directly in leaves
        if (vc->IsLeaf())
        {
                return vc->GetPvs().CountObjectsInPvs();
        }
        
        // 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->mPvsSize;
        }
        
        // 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();
        ObjectPvs newPvs;

        // sum different intersectables
        for (it = leaves.begin(); it != it_end; ++ it)
        {
                ObjectPvsMap::iterator oit, oit_end = (*it)->GetPvs().mEntries.end();

                for (oit = (*it)->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
                {
                        Intersectable *intersect = (*oit).first;

                        if (!intersect->Mailed())
                        {
                                intersect->Mail();
                                newPvs.AddSample(intersect, (*oit).second.mSumPdf);
                        }
                }
        }

        return newPvs.CountObjectsInPvs();
}


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)
        {
                ObjectPvsMap::iterator oit, oit_end = (*it)->GetPvs().mEntries.end();

                for (oit = (*it)->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
                {
                        Intersectable *intersect = (*oit).first;

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

        return pvsSize;
}


int ViewCellsTree::GetPvsSizeForCompressedStorage(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->mPvsSize;
        }

        // 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
                pvsSize += CountDiffPvs(root);
        }

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

        while (!tstack.empty())
        {
                vc = tstack.top();
                tstack.pop();
                // matt: bug! must evaluate kd pvss also
                pvsSize += CountDiffPvs(vc);

                if (!vc->IsLeaf())
                {
                        ViewCellInterior *interior = dynamic_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;
}


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 += CountDiffPvs(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 += CountDiffPvs(vc);

                if (!vc->IsLeaf())
                {
                        ViewCellInterior *interior = dynamic_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;
}


int ViewCellsTree::GetPvsSize(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 = GetPvsSizeForLeafStorage(vc);                 
                break;
        case COMPRESSED:
                pvsSize = GetPvsSizeForCompressedStorage(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().CountObjectsInPvs();     
                break;
        }

        return pvsSize;  
}


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
{
        const float entrySize = 
                sizeof(PvsData) + sizeof(Intersectable *);

        return (float)GetStoredPvsEntriesNum(vc) * entrySize;
}


int ViewCellsTree::GetStoredPvsEntriesNum(ViewCell *root) const
{
        int pvsSize = root->GetPvs().GetSize();

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

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

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

        return pvsSize;         
}


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


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


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

        // propagate pvs up
        while (viewCell->GetParent())
        {
                viewCell->GetParent()->GetPvs().Merge(vc->GetPvs());
                viewCell = viewCell->GetParent();
        }

        if (vc->IsLeaf())
                return;

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

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

                if (viewCell != vc)
                {
                        viewCell->GetPvs().Merge(vc->GetPvs());
                }

                if (!viewCell->IsLeaf())
                {
                        ViewCellInterior *interior = dynamic_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 = dynamic_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);
        
        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 = dynamic_cast<ViewCellInterior *>(vc);
                        costFunction.push_back(interior->GetMergeCost());

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

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

                }
        }
}


void ViewCellsTree::UpdateViewCellsStats(ViewCell *vc, 
                                                                                 ViewCellsStatistics &vcStat)
{
        ++ vcStat.viewCells;
                
        const int pvsSize = GetPvsSize(vc);

        vcStat.pvsSize += pvsSize;

        if (pvsSize == 0)
                ++ vcStat.emptyPvs;

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

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

        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 = dynamic_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)
{
        ObjectPvsMap::iterator it, it_end = viewCell->GetPvs().mEntries.end();

        for (it = viewCell->GetPvs().mEntries.begin(); it != it_end; ++ it)
        {
                stream << (*it).first->GetId() << " ";
        }
}


void ViewCellsTree::ExportViewCell(ViewCell *viewCell, 
                                                                   OUT_STREAM &stream, 
                                                                   const bool exportPvs)
{
        if (viewCell->IsLeaf())
        {
                stream << "<Leaf ";
                stream << "id=\"" << viewCell->GetId() << "\" ";
                stream << "active=\"" << dynamic_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 = dynamic_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 (0) 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 = dynamic_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";
}

}