source: GTP/trunk/Lib/Vis/Preprocessing/src/OspTree.cpp @ 1444

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

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


namespace GtpVisibilityPreprocessor {


#define USE_FIXEDPOINT_T 0


//-- static members

OspTree *OspTree::OspSubdivisionCandidate::sOspTree = NULL;

// variable for debugging volume contribution for heuristics
static float debugVol;

static bool ViewCellHasMultipleReferences(Intersectable *obj, 
                                                                                  ViewCell *vc, 
                                                                                  bool checkOnlyMailed)
{
        MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);
//return false;
        if (vdata)
        {
                // more than one view cell sees this object inside different kd cells
                if (!checkOnlyMailed || !vdata->Mailed())
                {
                        if (checkOnlyMailed)
                                vdata->Mail();
                        //Debug << "sumpdf: " << vdata->mSumPdf << endl;
                        if (vdata->mSumPdf > 1.5f) 
                                return true;
                }
        }

        return false;
}


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

        app << setprecision(4);

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

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

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

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

        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits[0] + splits[1] + splits[2] << endl;

        app << "#N_SPLITS ( Number of splits in axes x y z)\n";

        for (int i = 0; i < 3; ++ i)
                app << splits[i] << " ";
        app << endl;

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

        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n" 
                << minPvsNodes * 100 / (double)Leaves() << endl;

        //app << "#N_PMINOBJECTREFLEAVES  ( Percentage of leaves with minimal number of object ref)\n" 
                //<< minObjectNodes * 100 / (double)Leaves() << endl;

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

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

        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;

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

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

        app << "#N_INVALIDLEAVES (number of invalid leaves )\n" << invalidLeaves << endl;

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

//      app << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
        //app << "#N_PVS: " << pvs << endl;

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


/******************************************************************/
/*                  class VspNode implementation                  */
/******************************************************************/


/*******************************************************************/
/*                  class OspTree implementation                   */
/*******************************************************************/


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


OspTree::OspTree(const KdTree &kdTree)
{
        ReadEnvironment();

        // copy values from kd tree
        mCopyFromKdTree = true;
        mBoundingBox = kdTree.GetBox();
        mRoot = kdTree.GetRoot();

        mSubdivisionCandidates = new vector<SortableEntry>;
}


OspTree::~OspTree()
{
        // delete kd intersectables
        KdIntersectableMap::iterator it, it_end = mKdIntersectables.end();

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

        // delete hierarchy only if not using the hierarchy from 
        // some other kd tree (otherwise the other kd tree has to do the job)
        if (!mCopyFromKdTree)
                DEL_PTR(mRoot);

        DEL_PTR(mSubdivisionCandidates);
        mSubdivisionStats.close();
}


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

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

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

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

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

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


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

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

        //-- debug output

        Debug << "******* OSP tree options ******** " << endl;

    Debug << "max depth: " << mTermMaxDepth << endl;
        Debug << "min probabiliy: " << mTermMinProbability << endl;
        Debug << "min objects: " << mTermMinObjects << endl;
        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
        Debug << "miss tolerance: " << mTermMissTolerance << endl;
        Debug << "max leaves: " << mTermMaxLeaves << endl;

        Debug << "randomize: " << randomize << endl;

        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
        Debug << "max memory: " << mMaxMemory << endl;
        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
        
        Debug << "split borders: " << mSplitBorder << endl;
        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
        Debug << endl;
}


void OspTree::SplitObjects(KdLeaf *parent,
                                                   const AxisAlignedPlane& splitPlane,
                                                   const ObjectContainer &objects,
                                                   ObjectContainer &front,
                                                   ObjectContainer &back)
{
        ObjectContainer::const_iterator oit, oit_end = objects.end();

    for (oit = objects.begin(); oit != oit_end; ++ oit) 
        {
                Intersectable *object = *oit;
                
                // initialise leaf references of objects
                if (parent->IsRoot())
                {
                        object->mReferences = 1;
                }

                // remove parent ref
                -- object->mReferences;

                // determine the side of this ray with respect to the plane
                const AxisAlignedBox3 box = object->GetBox();

                if (box.Max(splitPlane.mAxis) >= splitPlane.mPosition)
                {
            front.push_back(object);
                        ++ object->mReferences;
                }

                if (box.Min(splitPlane.mAxis) < splitPlane.mPosition)
                {
                        back.push_back(object);
                        ++ object->mReferences;
                }
        }

        mOspStats.objectRefs -= (int)objects.size();
        mOspStats.objectRefs += (int)(back.size() + front.size());
}


KdInterior *OspTree::SubdivideNode(
                                                                   const AxisAlignedPlane &splitPlane,
                                                                   const OspTraversalData &tData,
                                                                   OspTraversalData &frontData,
                                                                   OspTraversalData &backData
                                                                   )
{
        KdLeaf *leaf = tData.mNode;
        
        // two new leaves
    mOspStats.nodes += 2;

        // add the new nodes to the tree
        KdInterior *node = new KdInterior(leaf->mParent);

        const int axis = splitPlane.mAxis;
        const float position = splitPlane.mPosition;

        node->mAxis = axis;
        node->mPosition = position;
        node->mBox = tData.mBoundingBox;


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

        // bounding boxes: split front and back node geometry
        tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition, 
                frontData.mBoundingBox, backData.mBoundingBox);

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

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

        //-- classify objects

        int objectsBack = 0;
        int objectsFront = 0;

    ObjectContainer::const_iterator mi, mi_end = leaf->mObjects.end();

        for ( mi = leaf->mObjects.begin(); mi != mi_end; ++ mi) 
        {
                // determine the side of this ray with respect to the plane
                const AxisAlignedBox3 box = (*mi)->GetBox();
                
                if (box.Max(axis) >= position) 
                        ++ objectsFront;
    
                if (box.Min(axis) < position)
                        ++ objectsBack;
        }


        // TODO matt: compute pvs
        //frontData.mPvs = objectsFront;
        //backData.mPvs = objectsBack;

        //CheckViewCellsPvs(leaf, *tData.mRays);
        RemoveParentViewCellsPvs(leaf, *tData.mRays);


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

        KdLeaf *back = new KdLeaf(node, objectsBack);
        KdLeaf *front = new KdLeaf(node, objectsFront);

        KdInterior *parent = leaf->mParent;

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

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

        SplitObjects(leaf, splitPlane, leaf->mObjects, front->mObjects, back->mObjects);
        
        FilterRays(front, *tData.mRays, *frontData.mRays);
        FilterRays(back, *tData.mRays, *backData.mRays);

        ++ mOspStats.splits[splitPlane.mAxis];

        //-- eval view cells pvs
        
        // remove parent pvs update pvs of left and right leaves
        // Note: It is necessary to update first left and then right pvs.
        // We need to store the view cells seen by each object,
        // but also if the view cells are seen as part of two different
        // kd leaves, which is stored in the pdf component of the pvs entry. 
        // Because if an object is part of a view cell more than once,
        // it cannot be removed from the pvs by splitting a kd node where
        // the view cell sees only the other child of the node.
        // This is important during the subdivision

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

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


        ProcessMultipleRefs(front);
    ProcessMultipleRefs(back);

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

        // compute probability, i.e., volume of seen view cells
        frontData.mProbability = EvalViewCellsVolume(front, *frontData.mRays);
        backData.mProbability =  EvalViewCellsVolume(back, *backData.mRays);

        //delete leaf;
        return node;
}


KdNode *OspTree::Subdivide(SplitQueue &tQueue,
                                                   SubdivisionCandidate *splitCandidate,
                                                   const bool globalCriteriaMet)
{
        OspSubdivisionCandidate *sc = dynamic_cast<OspSubdivisionCandidate *>(splitCandidate);
        OspTraversalData &tData = sc->mParentData;

        KdNode *newNode = tData.mNode;

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

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

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

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

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

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

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

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


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

                //-- store additional info
                if (mStoreRays)
                {
                        KdLeaf *leaf = dynamic_cast<KdLeaf *>(newNode);

                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();

                        for (it = tData.mRays->begin(); it != it_end; ++ it)
                        {
                                (*it).mRay->Ref();                      
                                leaf->mVssRays.push_back((*it).mRay);
                        }
                }
        }
        
        //-- cleanup
        tData.Clear();
        
        return newNode;
}


void OspTree::EvalSubdivisionCandidate(OspSubdivisionCandidate &splitCandidate)
{
        // compute locally best split plane
        const float ratio = SelectSplitPlane(splitCandidate.mParentData, 
                                                                                 splitCandidate.mSplitPlane);

        const bool success = ratio < mTermMaxCostRatio;
        float oldRenderCost;

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

        splitCandidate.SetRenderCostDecrease(renderCostDecr);

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

#endif

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

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


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


inline  bool OspTree::GlobalTerminationCriteriaMet(const OspTraversalData &data) const
{
        // matt: TODO
        return (0
                || (mOspStats.Leaves() >= mTermMaxLeaves) 
                //|| mOutOfMemory
                || (mGlobalCostMisses >= mTermGlobalCostMissTolerance)
                );
}


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

        /*if (data.mPvs > mOspStats.maxPvs)
        {
                mOspStats.maxPvs = data.mPvs;
        }

        mOspStats.pvs += data.mPvs;

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

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

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


float OspTree::EvalLocalCostHeuristics(const OspTraversalData &tData,
                                                                           const AxisAlignedBox3 &box,
                                                                           const int axis,
                                                                           float &position,
                                                                           int &objectsFront,
                                                                           int &objectsBack)
{
        RayInfoContainer usedRays;
        int mMaxTests = 500000; // HACK

        if (mMaxTests < (int)tData.mRays->size())
        {
                GetRayInfoSets(*tData.mRays, mMaxTests, usedRays);
        }
        else
        {
                usedRays = *tData.mRays;
        }

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

        const float sizeBox = maxBox - minBox;

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

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

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

        float totalVol = 0;
        float voll = 0;
        float volr = totalVol;

        ViewCellContainer touchedViewCells;
        
        const float totalRenderCost = PrepareHeuristics(tData, touchedViewCells);
        float renderCost = totalRenderCost;

        /// this is kind of a reverse pvs
        const int totalPvs = (int)tData.mNode->mObjects.size();
        
        int pvsl = 0;
        int pvsr = totalPvs;

        float sum = (float)totalVol * sizeBox;

        Debug << "total render cost: " << renderCost / viewSpaceVol << endl;

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

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

        float volBack = voll;
        float volFront = volr;

        int pvsBack = pvsl;
        int pvsFront = pvsr;
        
        float minRenderCost = 1e20f;

        debugVol = 0;



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

        
        //-- traverse through events and find best split plane
        
        vector<SortableEntry>::const_iterator ci, ci_end = mSubdivisionCandidates->end();
        //minRenderCost = RandomValue(0,viewSpaceVol);

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

                EvalHeuristicsContribution(tData.mNode, 
                                                                   *ci, 
                                                                   renderCost,
                                                                   touchedViewCells
                                                                  );

                // Note: sufficient to compare size of bounding boxes of front and back side?
                if (((*ci).mPos >= minBand) && ((*ci).mPos <= maxBand))
                {
                        float currentPos;
                        Debug << "current min: " << minRenderCost / viewSpaceVol << endl;

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

                        if (renderCost < minRenderCost)
                        {
                                splitPlaneFound = true;
                                minRenderCost = renderCost;
                                position = currentPos;
                                Debug << "pos: " << position << endl;
                        }
                }
        }
//splitPlaneFound = true;
        if (splitPlaneFound)
        {
                ratio = minRenderCost / totalRenderCost;
        }

        Debug << "\n§§§§ eval local cost §§§§" << endl
                  << "old rc: " << totalRenderCost / viewSpaceVol << " new rc: " << minRenderCost / viewSpaceVol << endl
                  << "render cost decrease: " << (totalRenderCost - minRenderCost) / viewSpaceVol  << endl;

        return ratio;
}


void OspTree::SortSubdivisionCandidates(const OspTraversalData &tData,
                                                                  const int axis, 
                                                                  float minBand, 
                                                                  float maxBand)

{
        mSubdivisionCandidates->clear();

        RayInfoContainer *rays = tData.mRays;
        KdLeaf *leaf = tData.mNode;

        int requestedSize = 2 * (int)rays->size() + 2 * (int)leaf->mObjects.size();

        // creates a sorted split candidates array
        if (mSubdivisionCandidates->capacity() > 500000 &&
                requestedSize < (int)(mSubdivisionCandidates->capacity() / 10) )
        {
        delete mSubdivisionCandidates;
                mSubdivisionCandidates = new vector<SortableEntry>;
        }

        mSubdivisionCandidates->reserve(requestedSize);

        float pos;

        //-- insert ray queries
        //-- we are intersested only in rays which intersect an object that
        //-- is part of the kd node because they can induce a change in view cell
        //-- volume on the left and the right part. 
        //-- Note that they cannot induce a change in pvs size!!

        RayInfoContainer::const_iterator rit, rit_end = rays->end();

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

                const bool positive = ray->HasPosDir(axis);
        
                if (EndPointInsideNode(leaf, *ray, true))
                {
                        pos = ray->mTermination[axis];

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

                // if hitpoint with object is inside this node
                if (EndPointInsideNode(leaf, *ray, false))
                {
                        pos = ray->mOrigin[axis];
        
                        mSubdivisionCandidates->push_back(
                                SortableEntry(SortableEntry::BOX_INTERSECT, 
                                                          pos, 
                                                          ray->mOriginObject, 
                                                          ray)
                                                          );
                }
        }


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

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

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

                mSubdivisionCandidates->push_back(
                        SortableEntry(SortableEntry::BOX_MIN,
                                                  box.Min(axis),
                                                  object,
                                                  NULL)
                                                  );
   
    
                mSubdivisionCandidates->push_back(
                        SortableEntry(SortableEntry::BOX_MAX,
                                                  box.Max(axis),
                                                  object,
                                                  NULL)
                                                  );
        }

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


const OspTreeStatistics &OspTree::GetStatistics() const
{
        return mOspStats;
}


void OspTree::PrepareHeuristics(const VssRay &ray, 
                                                                ViewCellContainer &touchedViewCells)
{
        // collect view cells and set mail + counter
        ViewCellContainer viewCells;
        mVspTree->GetViewCells(ray, viewCells);
        
        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

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

                if (!vc->Mailed())
                {
                        vc->Mail(); // set as belonging to front leaf
                        vc->mCounter = 0;
                        touchedViewCells.push_back(vc);
                }
                
                // view cell volume already added => just increase counter
                ++ vc->mCounter;
        }
}


float OspTree::PrepareHeuristics(const OspTraversalData &tData,
                                                                 ViewCellContainer &touchedViewCells)
{       
        float renderCost = 0;

        Intersectable::NewMail(3);
        ViewCell::NewMail(3);
        MailablePvsData::NewMail();

        KdLeaf *leaf = tData.mNode;

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

        for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit).mRay;
        PrepareHeuristics(*ray, touchedViewCells);
        }

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

        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *obj = *oit;
                obj->Mail(); // set as belonging to front leaf
                
                ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

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

                        if (!vdata) // should never happen
                                continue;

                        if (!vdata->Mailed())
                        {
                                vdata->Mail();
                                renderCost += vc->GetVolume();
                        }
                }
        }

        return renderCost;
}


void OspTree::AddObjectContribution(KdLeaf *leaf,
                                                  Intersectable *obj,
                                                  ViewCellContainer &touchedViewCells,
                                                  float &renderCost)
{
        //Debug << "add contri to obj: " << obj->mMailbox - Intersectable::sMailId << endl;
        obj->Mail(2); // set as belonging to both leafs

        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();
        
        for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = *vit;
        
                // if obj not previously associated with this view cell => increase render cost
                if (vc->Mailed(1) && !ViewCellHasMultipleReferences(obj, vc, false))
                {       
                        //Debug << "add to rendercost: " << vc->GetVolume() / mVspTree->GetBoundingBox().GetVolume() << endl;
                        renderCost += vc->GetVolume();
                }
        }
}


void OspTree::SubtractObjectContribution(KdLeaf *leaf,
                                                                Intersectable * obj,
                                                                ViewCellContainer &touchedViewCells,
                                                                float &renderCost)
{
        //Debug << "subtract contri from obj: " << obj->mMailbox - Intersectable::sMailId << endl;
    obj->Mail(1); // set as belonging to back leaf
        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();
        
        for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = *vit;

                // if obj was previously associated with this view cell but is not now 
                // => decrease render cost
                if (vc->Mailed() &&     !ViewCellHasMultipleReferences(obj, vc, false))
                {
                        //Debug << "subtract from rendercost: " << vc->GetVolume() / mVspTree->GetBoundingBox().GetVolume() << endl;
                        renderCost -= vc->GetVolume();
                }
        }
}


void OspTree::EvalRayContribution(KdLeaf *leaf,
                                                                  const VssRay &ray, 
                                                                  float &renderCost)
{
        ViewCellContainer viewCells;

        mVspTree->GetViewCells(ray, viewCells);

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


void OspTree::EvalViewCellContribution(KdLeaf *leaf, 
                                                                           ViewCell *viewCell,
                                                                           float &renderCost)
{
        //Debug << "**************" << endl;
        //Debug << "vc contri: " << renderCost / mVspTree->GetBoundingBox().GetVolume() << " " << viewCell->mMailbox - ViewCell::sMailId << " counter: " << viewCell->mCounter << endl;

        const float vol = viewCell->GetVolume();

        if (viewCell->Mailed())
        {
                viewCell->Mail(2); // view cell can be seen from both nodes

        // we now see view cell from both nodes => add contri
                ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

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

                        // was render cost already added?
                        if (!ViewCellHasMultipleReferences(obj, viewCell, false) &&
                                obj->Mailed(1))
                        {
                                //Debug << "obj mail 1!! " << vol / mVspTree->GetBoundingBox().GetVolume() << endl;
                                renderCost += vol;
                        }
                }
        }

        if (-- viewCell->mCounter == 0)
        {
                viewCell->Mail(1); // view cell can be seen from back node only

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

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

                        // can render cost be be reduced?
                        if (!ViewCellHasMultipleReferences(obj, viewCell, false) &&     obj->Mailed())
                        {
                                renderCost -= vol;
                        }
                }
        }
        //Debug << "new rc: " << renderCost / mVspTree->GetBoundingBox().GetVolume() << " counter2: " << viewCell->mCounter << endl;
        //Debug << "**************"<<endl;
}


void OspTree::EvalHeuristicsContribution(KdLeaf *leaf,
                                                                                 const SortableEntry &ci,
                                                                                 float &renderCost,
                                                                                 ViewCellContainer &touchedViewCells)
{
        Intersectable *obj = ci.mObject;
        VssRay *ray = ci.mRay;

        // switch between different types of events
        switch (ci.mType) 
        {
                case SortableEntry::BOX_MIN:
                        cout << "<";
                        AddObjectContribution(leaf, ci.mObject, touchedViewCells, renderCost);  
                        break;
                        
                case SortableEntry::BOX_MAX:
                        cout << ">";
                        SubtractObjectContribution(leaf, ci.mObject, touchedViewCells, renderCost);
                        break;

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

                default:
                        cout << "should not come here" << endl;
                        break;
        }

        //cout << "vol left: " << volLeft << " vol right " << volRight << endl;
}


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


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


float OspTree::SelectSplitPlane(const OspTraversalData &tData,
                                                                AxisAlignedPlane &plane)
{
        float nPosition[3];
        float nCostRatio[3];

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

        if (mOnlyDrivingAxis)
        {
                sAxis = box.Size().DrivingAxis();
        }
        

        // -- evaluate split cost for all three axis
        for (int axis = 0; axis < 3; ++ axis)
        {
                if (!mOnlyDrivingAxis || (axis == sAxis))
                {
                        if (1 || mUseCostHeuristics)
                        {
                                //-- place split plane using heuristics
                                int objectsFront, objectsBack;

                                nCostRatio[axis] =
                                        EvalLocalCostHeuristics(tData,
                                                                           tData.mBoundingBox,
                                                                           axis,
                                                                           nPosition[axis],
                                                                           objectsFront,
                                                                           objectsBack);                        
                        }
                        /*      else 
                        {
                                //-- split plane position is spatial median

                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;

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

        //-- assign values

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

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


bool OspTree::EndPointInsideNode(KdLeaf *leaf, 
                                                                 VssRay &ray, 
                                                                 bool isTermination) const
{
        // test if the leaf where the hitpoint is located is the current leaf
        if (isTermination)
        {
                 return ray.mTerminationObject && (GetLeaf(ray.mTermination, ray.mTerminationNode) == leaf);
        }
        else
        {
                return false; // no origin object!
                return ray.mOriginObject && (GetLeaf(ray.mOrigin, ray.mOriginNode) == leaf);
        }
}


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

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


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


int OspTree::ClassifyRay(VssRay *ray, 
                                                 KdLeaf *leaf, 
                                                 const AxisAlignedPlane &plane) const
{
        const bool originInside = EndPointInsideNode(leaf, *ray, false);
    const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

        const bool originGt = (ray->mOrigin[plane.mAxis] >= plane.mPosition);
        const bool originLt = (ray->mOrigin[plane.mAxis] <= plane.mPosition);

        const bool terminationGt = (ray->mTermination[plane.mAxis] >= plane.mPosition);
        const bool terminationLt = (ray->mTermination[plane.mAxis] <= plane.mPosition);

        // add view cell volume to front volume
        const bool addToFront = 
                ((originInside && originGt) || (terminationInside && terminationGt));

        // add view cell volume to back volume
        const bool addToBack = 
                ((originInside && originLt/*!originGt*/) || (terminationInside && terminationLt/*terminationGt*/));

        // classify ray with respect to the child nodes the view cells contribute
        if (addToFront && addToBack)
                return 0;
        else if (addToBack)
                return -1;
        else if (addToFront)
                return 1;

        return -2;
}


int OspTree::ClassifyRays(const RayInfoContainer &rays, 
                                                  KdLeaf *leaf, 
                                                  const AxisAlignedPlane &plane,
                                                  RayInfoContainer &frontRays,
                                                  RayInfoContainer &backRays) const
{
        RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

                bool originGt = false;
                bool terminationGt = false;
                
                Intersectable *obj = ray->mOriginObject;

                if (0 && obj)
                {
                        originGt = ray->mOrigin[plane.mAxis] >= plane.mPosition;
                }

                obj = ray->mTerminationObject;
                
                if (obj)
                {
                        terminationGt = ray->mTermination[plane.mAxis] >= plane.mPosition;
                }

                // either begin or end point on front side
                const bool addToFront = originGt || terminationGt;

                // add view cell volume to back volume
                const bool addToBack = !originGt || !terminationGt;
        
                // classify ray with respect to the child nodes the view cells contribute
                if (addToFront)
                        frontRays.push_back(*rit);
                if (addToBack)
                        backRays.push_back(*rit);
        }

        return 0;
}


float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
                                                                          const OspTraversalData &tData,
                                                                          float &normalizedOldRenderCost) const
{
        float pvsFront = 0;
        float pvsBack = 0;

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

        Intersectable::NewMail(3);
        KdLeaf *leaf = tData.mNode;

        const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();
        const int totalPvs = (int)leaf->mObjects.size();
        
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

        // evaluate reverse pvs and view cell volume on left and right cell
        // note: should I take all leaf objects or rather the objects hit by rays?
        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                int pvsContri = 0;

                Intersectable *obj = *oit;
                const AxisAlignedBox3 box = obj->GetBox();

                // test if box falls in left / right child node
                if (box.Max(candidatePlane.mAxis) >= candidatePlane.mPosition)
                {
                        ++ pvsFront;
                        obj->Mail();            
                }

                if (box.Min(candidatePlane.mAxis) < candidatePlane.mPosition)
                {
                        ++ pvsBack;

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

        ViewCellContainer touchedViewCells;

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

        RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
        //map<Intersectable *, ViewCellContainer> objectsMap;

        for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
        {
                VssRay *ray = (*rit).mRay;

                // add volume to volumes of left and / or right children
                // if one of the ray end points is inside
                const int classification = ClassifyRay(ray, leaf, candidatePlane);

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

                // traverse through view cells and classify them according
                // to them being seen from to back / front / front and back node
                for (vit = viewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;

                        // if not previously mailed
                        if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
                        {
                                touchedViewCells.push_back(vc);
                        }

                        // classify / mail the view cell
                        MailViewCell(vc, classification);
                }
        }

        // evaluate view cells volume contribution 
        // with respect to the mail box classification
        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

        for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
        {
                AddViewCellVolumeContri(*vit, pFront, pBack, pFrontAndBack);
        }

        //////////////////////////////////////////////
        //
        // evaluate contribution of objects which are part of other kd nodes
        // which are seen by the same view cell
        // These contributions cannot be reduced by splitting, because
        // the object would still be part of the pvs

        float newRc = 0;
        float rc = 0;

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

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

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

                        if (!vdata)
                        {
                                Debug << "error!!" << endl;
                                continue;
                        }

                        if (!vdata->Mailed())
                        {
                                vdata->Mail();
                                rc += vc->GetVolume();

                                if ((vdata->mSumPdf > 1.5) || 
                                        vc->Mailed(2) || 
                                        obj->Mailed(2) ||
                                        (obj->Mailed() && vc->Mailed()) ||
                                        (obj->Mailed(1) && vc->Mailed(1))
                                        )
                                {
                                        newRc += vc->GetVolume();
                                }
                        }
                }
        }


        /////////////////////////////
        //-- pvs rendering heuristics

        const float oldRenderCost = rc;
        
        // sum up the terms: 
        // the view cells seeing
        // a) the left node are weighted by the #left node objects
        // b) the right node are weighted by the #right node objects
        // c) both nodes are weighted by the #parent node objects
        const float newRenderCost = newRc;


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

        Debug << "\nrender cost decrease: " << endl
                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
                  << "render cost decrease: " << renderCostDecrease << endl;

        normalizedOldRenderCost = oldRenderCost / viewSpaceVol;

                
        return renderCostDecrease;
}


void OspTree::ComputeBoundingBox(const ObjectContainer &objects) 
{
        mBoundingBox.Initialize();

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

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

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


void OspTree::ProcessMultipleRefs(KdLeaf *leaf) const
{
        // find objects from multiple kd-leaves
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
        {
                Intersectable *object = *oit;
                
                if (object->mReferences > 1)
                {
                        leaf->mMultipleObjects.push_back(object);
                }
        }
}


void OspTree::CollectLeaves(vector<KdLeaf *> &leaves) const
{
        stack<KdNode *> nodeStack;
        nodeStack.push(mRoot);

        while (!nodeStack.empty()) 
        {
                KdNode *node = nodeStack.top();
                nodeStack.pop();
                if (node->IsLeaf()) 
                {
                        KdLeaf *leaf = (KdLeaf *)node;
                        leaves.push_back(leaf);
                }
                else 
                {
                        KdInterior *interior = (KdInterior *)node;
                        nodeStack.push(interior->mBack);
                        nodeStack.push(interior->mFront);
                }
        }
}


AxisAlignedBox3 OspTree::GetBoundingBox(KdNode *node) const
{
        if (!node->mParent)
                return mBoundingBox;

        if (!node->IsLeaf())
        {
                return (dynamic_cast<KdInterior *>(node))->mBox;                
        }

        KdInterior *parent = dynamic_cast<KdInterior *>(node->mParent);

        AxisAlignedBox3 box(parent->mBox);

        if (parent->mFront == node)
                box.SetMin(parent->mAxis, parent->mPosition);
    else
                box.SetMax(parent->mAxis, parent->mPosition);

        return box;
}


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

        ViewCell::NewMail();

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

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

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

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


void OspTree::CollectDirtyCandidates(OspSubdivisionCandidate *sc, 
                                                                         vector<SubdivisionCandidate *> &dirtyList)
{
        OspTraversalData &tData = sc->mParentData;
        KdLeaf *node = tData.mNode;
        
        ViewCell::NewMail();
        ViewCellContainer viewCells;
        ViewCellContainer tmpViewCells;

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

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

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

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

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


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

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

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


KdNode *OspTree::GetRoot() const
{
        return mRoot;
}


KdLeaf *OspTree::GetLeaf(const Vector3 &pt, KdNode *node) const
{
        // start from root of tree
        if (node == NULL)
        {
                node = mRoot;
        }

        stack<KdNode *> nodeStack;
        nodeStack.push(node);
  
        KdLeaf *leaf = NULL;
  
        while (!nodeStack.empty())  
        {
                KdNode *node = nodeStack.top();
                nodeStack.pop();
        
                if (node->IsLeaf()) 
                {
                        leaf = dynamic_cast<KdLeaf *>(node);
                } 
                else    
                {       
                        // find point
                        KdInterior *interior = dynamic_cast<KdInterior *>(node);
        
                        if (interior->mPosition > pt[interior->mAxis])
                        {
                                nodeStack.push(interior->mBack);
                        }
                        else
                        {
                                nodeStack.push(interior->mFront);
                        }
                }
        }
  
        return leaf;
}


void OspTree::PreprocessRays(KdLeaf *root,
                                                         const VssRayContainer &sampleRays,
                                                         RayInfoContainer &rays)
{
        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
        RayInfoContainer clippedRays;

        long startTime = GetTime();

        cout << "storing rays ... ";

        Intersectable::NewMail();

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

                float minT, maxT;
                static Ray hray;

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

                        clippedRays.push_back(RayInfo(ray, minT / len, maxT / len));

                        // HACK: reset nodes for the case we have used object kd tree for
                        // tree building.
                        ray->mOriginNode = NULL;
                        ray->mTerminationNode = NULL;
                }
        }

        FilterRays(root, clippedRays, rays);

        cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
}



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

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

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

                // get classification and receive new t
                //-- test if start point behind or in front of plane
                const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
                        
                if (side == 0)
                {
                        ++ splits;

                        if (ray->HasPosDir(plane.mAxis))
                        {
                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
                        }
                        else
                        {
                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
                        }
                }
                else if (side == 1)
                {
                        frontRays.push_back(bRay);
                }
                else
                {
                        backRays.push_back(bRay);
                }
        }

        return splits;
}


int OspTree::FilterRays(KdLeaf *leaf, 
                                                const RayInfoContainer &rays,
                                                RayInfoContainer &filteredRays)
{
        RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

                // test if intersection point with one of the objects is inside this node
                const bool originInside = EndPointInsideNode(leaf, *ray, false);
        const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

                if (originInside || terminationInside)
                {
                        filteredRays.push_back(ray);
                }
        }

        return 0;
}
                                         

int OspTree::CheckViewCellsPvs(KdLeaf *leaf, 
                                                           const RayInfoContainer &rays) const
{
        RayInfoContainer::const_iterator rit, rit_end = rays.end();

        Intersectable::NewMail();
        ObjectContainer touchedObjects;
        

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

                if (ray->mTerminationObject)
                {
                        Intersectable *obj = ray->mTerminationObject;

                        if (!obj->Mailed())
                        {
                                obj->Mail();
                                touchedObjects.push_back(obj);
                        }
                }

                if (ray->mOriginObject)
                {
                        Intersectable *obj = ray->mOriginObject;

                        if (!obj->Mailed())
                        {
                                obj->Mail();
                                touchedObjects.push_back(obj);
                        }
                }
        }
        Debug << "touched view cells: " << (int)touchedObjects.size() << endl;
        ObjectContainer::const_iterator it, it_end = touchedObjects.end();
        for (it = touchedObjects.begin(); it != it_end; ++ it)
        {
                Debug << "\nobj: " << (*it) << " vc pvs size: " << (*it)->mViewCellPvs.GetSize() << endl << endl;
                ViewCellPvsMap::const_iterator mit, mit_end = (*it)->mViewCellPvs.mEntries.end();

                for (mit = (*it)->mViewCellPvs.mEntries.begin(); mit != mit_end; ++ mit)
                {
                        Debug << "new sumpdf: " << (*mit).second.mSumPdf << endl;
                }
        }

        return 0;
}


KdIntersectable *OspTree::GetOrCreateKdIntersectable(KdNode *node)
{
        // search nodes
        std::map<KdNode *, KdIntersectable *>::
                const_iterator it = mKdIntersectables.find(node);

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

        // not in map => create new entry
        KdIntersectable *kdObj = new KdIntersectable(node);
        mKdIntersectables[node] = kdObj;

        return kdObj;
}


/*
void OspTree::AddViewCellVolume(ViewCell *vc,
                                                                const int cf,
                                                                float &frontVol,
                                                                float &backVol,
                                                                float &totalVol) const
{
        //const float renderCost = mViewCellsManager->SimpleRay &raynderCost(obj);
        const float vol = vc->GetVolume();

        // view cell not found yet => new
        if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
        {
                totalVol += vol;
        }

        if (cf >= 0) // front volume
        {
                if (!vc->Mailed() && !vc->Mailed(2))
                {
                        frontVol += vol;
                
                        // already in back volume => in both volumes
                        if (vc->Mailed(1))
                                vc->Mail(2);
                        else
                                vc->Mail();
                }
        }

        if (cf <= 0) // back volume
        {
                if (!vc->Mailed(1) && !vc->Mailed(2))
                {
                        backVol += vol;
                
                        // already in front volume => in both volume
                        if (vc->Mailed())
                                vc->Mail(2);
                        else
                                vc->Mail(1);
                }
        }
}
*/


void OspTree::MailViewCell(ViewCell *vc, const int cf) const
{
        if (vc->Mailed(2)) // already classified
                return;

        if (cf >= 0) // front volume
        {
                // already in back volume => in both volumes
                if (vc->Mailed(1))
                {
                        vc->Mail(2);
                }
                else
                {
                        vc->Mail();
                }
        }

        if (cf <= 0) // back volume
        {
                // already in front volume => in both volume
                if (vc->Mailed())
                {
                        vc->Mail(2);
                }
                else
                {
                        vc->Mail(1);
                }
        }
}


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


float OspTree::EvalViewCellsVolume(KdLeaf *leaf, 
                                                                   const RayInfoContainer &rays) const
{
        float vol = 0;
        ViewCell::NewMail();

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

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

                ViewCellContainer viewCells;
                mVspTree->GetViewCells(*ray, viewCells);

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

                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                vol += vc->GetVolume();
                        }
                }               
        }

        return vol;
}


bool OspTree::AddViewCellToObjectPvs(Intersectable *obj, 
                                                                         ViewCell *vc, 
                                                                         float &contribution,
                                                                         bool onlyMailed) const
{       
        contribution = 0; // todo

        MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

        if (!vdata)
        {
                vdata = obj->mViewCellPvs.AddSample2(vc, 1);
        }
        else if (!onlyMailed || !vdata->Mailed())
        {
                obj->mViewCellPvs.AddSample(vc, 1);
        }
        
        vdata->Mail();

        return true;
}


void OspTree::CollectTouchedViewCells(const RayInfoContainer &rays, 
                                                                          ViewCellContainer &touchedViewCells) const
{
        ViewCell::NewMail();
        
        RayInfoContainer::const_iterator rit, rit_end = rays.end();

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

                ViewCellContainer viewCells;
                mVspTree->GetViewCells(*ray, viewCells);

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

                // traverse through view cells and classify them according
                // to them being seen from to back / front / front and back node
                for (vit = viewCells.begin(); vit != vit_end; ++ vit)
                {
                        ViewCell *vc = *vit;
                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                touchedViewCells.push_back(vc);
                        }
                }
        }
}


int OspTree::UpdateViewCellsPvs(KdLeaf *leaf, 
                                                                const RayInfoContainer &rays) const

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

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

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

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

        return 0;
}


int OspTree::RemoveParentViewCellsPvs(KdLeaf *leaf, 
                                                                          const RayInfoContainer &rays
                                                                          ) const

{
        MailablePvsData::NewMail();

        ViewCellContainer touchedViewCells;
        CollectTouchedViewCells(rays, touchedViewCells);

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

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

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

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

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

        return 0;
}


float OspTree::EvalRenderCost(const VssRayContainer &myrays)
{
        float rc = 0;
        float prop = mVspTree->GetBoundingBox().GetVolume();    

        KdLeafContainer leaves;
        CollectLeaves(leaves);

        ViewCellContainer vcleaves;
        mVspTree->CollectViewCells(vcleaves, false);
        

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

        for (vit = vcleaves.begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = *vit;
                vc->GetPvs().Clear();
        }

        KdLeafContainer::const_iterator kit, kit_end = leaves.end();

        MailablePvsData::NewMail();

        for (kit = leaves.begin(); kit != kit_end; ++ kit)
        {
                KdLeaf *leaf = *kit;
                float newCost = 0;
                float vol = 0;

                ViewCell::NewMail();
                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
                ViewCellContainer touchedViewCells;

                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
                {
                        VssRay *ray = *rit;
                        
                        // test if intersection point with one of the objects is inside this node
                        const bool originInside = EndPointInsideNode(leaf, *ray, false);
                        const bool terminationInside = EndPointInsideNode(leaf, *ray, true);

                        if (originInside || terminationInside)
                        {
                                ViewCellContainer viewCells;
                                mVspTree->GetViewCells(*ray, viewCells);
                        
                                ViewCellContainer::const_iterator vit, vit_end = viewCells.end();

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

                                        // if not previously mailed
                                        if (!vc->Mailed())
                                        {
                                                vc->Mail();
                                                touchedViewCells.push_back(vc);
                                        }

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

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

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

                                if (!vdata)
                                {
                                        vc->GetPvs().AddSample(obj, 1);
                                        newCost += vc->GetVolume();
                                }

/*                              MailablePvsData *vdata = obj->mViewCellPvs.Find(vc);

                                if (!vdata || !vdata->Mailed())
                                {
                                        newCost += vc->GetVolume();
                                        if (!vdata)
                                                vdata = obj->mViewCellPvs.AddSample2(vc, 1);
                                        vdata->Mail();
                                }*/
                        }
                }

                rc += newCost;
        }

        return rc / prop;
}


float OspTree::EvalLeafCost(const OspTraversalData &tData)
{
        KdLeaf *leaf = tData.mNode;

        float rc = 0;
        //float prop = mVspTree->GetBoundingBox().GetVolume();  
        //float vol = 0;
        //ViewCell::NewMail();
        
        RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
        ViewCellContainer touchedViewCells;

        for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
        {
                VssRay *ray = (*rit).mRay;
                        
                // test if intersection point with one of the objects is inside this node

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

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

                        // if not previously mailed
                        if (!vc->Mailed())
                        {
                                vc->Mail();
                                touchedViewCells.push_back(vc);
                        }
                }
        }

        // clear pvs of involved view cells
        ViewCellContainer::const_iterator vit, vit_end = touchedViewCells.end();

        for (vit = touchedViewCells.begin(); vit != vit_end; ++ vit)
        {
                ViewCell *vc = *vit;
                vc->GetPvs().Clear();
        }
        
        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();

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

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

                        if (!vdata)
                        {
                                vc->GetPvs().AddSample(obj, 1);
                                rc += vc->GetVolume();
                                //prop += vc->GetVolume();
                        }
                }
        }
        
        return rc;
}


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

        return true;
}


void OspTree::ExportNode(KdNode *node, OUT_STREAM &stream)
{
        if (node->IsLeaf())
        {
                KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);

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

                ExportNode(interior->mBack, stream);
                ExportNode(interior->mFront, stream);

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


struct KdObjectsTraversalData
{
        KdNode *node;
        ObjectContainer *objects;
};


void OspTree::InsertObjects(KdNode *node, const ObjectContainer &objects)
{
        stack<KdObjectsTraversalData> tStack;

        while (!tStack.empty())
        {
                KdObjectsTraversalData tData = tStack.top();
        tStack.pop();

                KdNode *node = tData.node;
                
                if (node->IsLeaf())
                {
                        KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);

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

                        for (oit = tData.objects->begin(); oit != oit_end; ++ oit)
                        {
                                leaf->mObjects.push_back(*oit);
                        }
                }
                else // interior
                {
                        KdObjectsTraversalData frontData, backData;
                        KdInterior *interior = dynamic_cast<KdInterior *>(node);

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

                        ObjectContainer::const_iterator oit, oit_end = tData.objects->end();
                        
                    for (oit = tData.objects->begin(); oit != oit_end; ++ oit) 
                        {
                                Intersectable *object = *oit;
                
                                // determine the side of this ray with respect to the plane
                                const AxisAlignedBox3 box = object->GetBox();

                                if (box.Max(interior->mAxis) >= interior->mPosition)
                                {
                                        frontData.objects->push_back(object);
                                }

                                if (box.Min(interior->mAxis) < interior->mPosition)
                                {
                                        backData.objects->push_back(object);
                                }
                        }

                        tStack.push(backData);
                        tStack.push(frontData);
                }

                DEL_PTR(tData.objects);
        }
}


SubdivisionCandidate * OspTree::PrepareConstruction(const VssRayContainer &sampleRays,
                                                                                                        const ObjectContainer &objects,
                                                                                                        RayInfoContainer &rays)
{
        // store pointer to this tree
        OspSubdivisionCandidate::sOspTree = this;
        mOspStats.nodes = 1;
        
        // compute bounding box from objects
        ComputeBoundingBox(objects);

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

        //-- create new root
        KdLeaf *kdleaf = new KdLeaf(NULL, 0);
        kdleaf->mObjects = objects;
        mRoot = kdleaf;
        
        // get clipped rays
        PreprocessRays(kdleaf, sampleRays, rays);

        // probabilty is voume of all "seen" view cells
#if 1
        const float prop = EvalViewCellsVolume(kdleaf, rays);
#else
        const float prop = GetBoundingBox().GetVolume();
#endif

        //-- add first candidate for object space partition

        // create osp traversal data
        OspTraversalData oData(kdleaf, 0, &rays, 0, prop, mBoundingBox);

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

        UpdateViewCellsPvs(kdleaf, rays);

        EvalSubdivisionCandidate(*oSubdivisionCandidate);

        mTotalCost = (float)objects.size() * prop / 
                mVspTree->GetBoundingBox().GetVolume();

        EvalSubdivisionStats(*oSubdivisionCandidate);

        return oSubdivisionCandidate;
}


bool OspTree::AddLeafToPvs(KdLeaf *leaf, 
                                                   ViewCell *vc, 
                                                   const float pdf, 
                                                   float &contribution)
{
        // add kd intersecable to pvs
        KdIntersectable *kdObj = GetOrCreateKdIntersectable(leaf);
        
        return vc->AddPvsSample(kdObj, pdf, contribution);
}


}