source: trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp @ 475

#include "Plane3.h"
#include "VspBspTree.h"
#include "Mesh.h"
#include "common.h"
#include "ViewCell.h"
#include "Environment.h"
#include "Polygon3.h"
#include "Ray.h"
#include "AxisAlignedBox3.h"
#include <stack>
#include <time.h>
#include <iomanip>
#include "Exporter.h"
#include "Plane3.h"
#include "ViewCellBsp.h"

//-- static members
/** Evaluates split plane classification with respect to the plane's 
        contribution for a minimum number of ray splits.
*/
const float VspBspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
/** Evaluates split plane classification with respect to the plane's 
        contribution for balanced rays.
*/
const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};


int VspBspTree::sFrontId = 0; 
int VspBspTree::sBackId = 0;
int VspBspTree::sFrontAndBackId = 0;


/****************************************************************/
/*                  class VspBspTree implementation             */
/****************************************************************/

VspBspTree::VspBspTree(): 
mRoot(NULL),
mPvsUseArea(true),
mCostNormalizer(Limits::Small)
{
        mRootCell = new BspViewCell();

        Randomize(); // initialise random generator for heuristics

        //-- termination criteria for autopartition
        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
        environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);     
        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
        environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);

        //-- factors for bsp tree split plane heuristics
        environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
        environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
        environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);

        //-- termination criteria for axis aligned split
        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAxisAlignedCtDivCi);
        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
        
        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays", 
                                                         mTermMinRaysForAxisAligned);
        
        //-- partition criteria
        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
        environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mAxisAlignedSplitBorder);

        environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
        environment->GetIntValue("VspBspTree.maxTests", mMaxTests);

        // maximum number of view cells
        environment->GetIntValue("ViewCells.maxViewCells", mMaxViewCells);

        //--
        Debug << "******* VSP BSP options ******** " << endl;
    Debug << "max depth: " << mTermMaxDepth << endl;
        Debug << "min PVS: " << mTermMinPvs << endl;
        Debug << "min area: " << mTermMinArea << endl;
        Debug << "min rays: " << mTermMinRays << endl;
        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
        //Debug << "VSP BSP mininam accumulated ray lenght: ", mTermMinAccRayLength) << endl;
        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
        Debug << "miss tolerance: " << mTermMissTolerance << endl;
        Debug << "max view cells: " << mMaxViewCells << endl;
        Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
        Debug << "max plane candidates: " << mMaxRayCandidates << endl;
        
        Debug << "Split plane strategy: ";
        if (mSplitPlaneStrategy & RANDOM_POLYGON)
        {
                Debug << "random polygon ";
        }
        if (mSplitPlaneStrategy & AXIS_ALIGNED)
        {
                Debug << "axis aligned ";
        }
        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
        {
                mCostNormalizer += mLeastRaySplitsFactor;
                Debug << "least ray splits ";
        }
        if (mSplitPlaneStrategy & BALANCED_RAYS)
        {
                mCostNormalizer += mBalancedRaysFactor;
                Debug << "balanced rays ";
        }
        if (mSplitPlaneStrategy & PVS)
        {
                mCostNormalizer += mPvsFactor;
                Debug << "pvs";
        }

        Debug << endl;
}


const BspTreeStatistics &VspBspTree::GetStatistics() const 
{
        return mStat;
}


VspBspTree::~VspBspTree()
{
        DEL_PTR(mRoot);
        DEL_PTR(mRootCell);
}

int VspBspTree::AddMeshToPolygons(Mesh *mesh, 
                                                                  PolygonContainer &polys, 
                                                                  MeshInstance *parent)
{
        FaceContainer::const_iterator fi;
        
        // copy the face data to polygons
        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
        {
                Polygon3 *poly = new Polygon3((*fi), mesh);
                
                if (poly->Valid(mEpsilon))
                {
                        poly->mParent = parent; // set parent intersectable
                        polys.push_back(poly);
                }
                else
                        DEL_PTR(poly);
        }
        return (int)mesh->mFaces.size();
}

int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells, 
                                                          PolygonContainer &polys, 
                                                          int maxObjects)
{
        int limit = (maxObjects > 0) ? 
                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
  
        int polysSize = 0;

        for (int i = 0; i < limit; ++ i)
        {
                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
                {
                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
                        polysSize += 
                                AddMeshToPolygons(viewCells[i]->GetMesh(), 
                                                                  polys, 
                                                                  viewCells[i]);
                }
        }

        return polysSize;
}

int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects, 
                                                                 PolygonContainer &polys, 
                                                                 int maxObjects)
{
        int limit = (maxObjects > 0) ? 
                Min((int)objects.size(), maxObjects) : (int)objects.size();
  
        for (int i = 0; i < limit; ++i)
        {
                Intersectable *object = objects[i];//*it;
                Mesh *mesh = NULL;

                switch (object->Type()) // extract the meshes
                {
                case Intersectable::MESH_INSTANCE:
                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
                        break;
                case Intersectable::VIEW_CELL:
                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
                        break;
                        // TODO: handle transformed mesh instances
                default:
                        Debug << "intersectable type not supported" << endl;
                        break;
                }
                
        if (mesh) // copy the mesh data to polygons
                {
                        mBox.Include(object->GetBox()); // add to BSP tree aabb
                        AddMeshToPolygons(mesh, polys, mRootCell);
                }
        }

        return (int)polys.size();
}

void VspBspTree::Construct(const VssRayContainer &sampleRays)
{
    mStat.nodes = 1;
        mBox.Initialize();      // initialise BSP tree bounding box
        
        PolygonContainer polys;
        RayInfoContainer *rays = new RayInfoContainer();

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

        long startTime = GetTime();

        cout << "Extracting polygons from rays ...";

        Intersectable::NewMail();

        //-- extract polygons intersected by the rays
        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
        {
                VssRay *ray = *rit;
        
                if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
                {
                        ray->mTerminationObject->Mail();
                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
                }

                if (ray->mOriginObject && !ray->mOriginObject->Mailed())
                {
                        ray->mOriginObject->Mail();
                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
                }
        }

        // compute bounding box
        Polygon3::IncludeInBox(polys, mBox);

        //-- store rays
        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
        {
                VssRay *ray = *rit;
                
                float minT, maxT;

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

        mTermMinArea *= mBox.SurfaceArea(); // normalize
        mStat.polys = (int)polys.size();

        cout << "finished" << endl;

        Debug << "\nPolygon extraction: " << (int)polys.size() << " polys extracted from " 
                  << (int)sampleRays.size() << " rays in "
                  << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;

        Construct(polys, rays);

        // clean up polygons
        CLEAR_CONTAINER(polys);
}

void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
{
        VspBspTraversalStack tStack;

        mRoot = new BspLeaf();

        // constrruct root node geometry
        BspNodeGeometry *geom = new BspNodeGeometry();
        ConstructGeometry(mRoot, *geom);

        VspBspTraversalData tData(mRoot, 
                                                          new PolygonContainer(polys), 
                                                          0,
                                                          rays, 
                              ComputePvsSize(*rays), 
                                                          geom->GetArea(), 
                                                          geom);

        tStack.push(tData);

        mStat.Start();
        cout << "Contructing vsp bsp tree ... ";

        long startTime = GetTime();
        
        while (!tStack.empty()) 
        {
                tData = tStack.top();

            tStack.pop();
                // subdivide leaf node
                BspNode *r = Subdivide(tStack, tData);

                if (r == mRoot)
                        Debug << "VSP BSP tree construction time spent at root: " 
                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl << endl;
        }

        cout << "finished\n";

        mStat.Stop();
}

bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data, 
                                                                                const int numLeaves) const
{
        return 
                (((int)data.mRays->size() <= mTermMinRays) ||
                 (data.mPvs <= mTermMinPvs)   ||
                 (data.mArea <= mTermMinArea) ||
                 ((mStat.nodes / 2 + 1) >= mMaxViewCells) ||
                // (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
                 (data.mDepth >= mTermMaxDepth));
}

BspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, 
                                                           VspBspTraversalData &tData)
{
        BspNode *newNode = tData.mNode;

        if (!TerminationCriteriaMet(tData, (int)tStack.size()))
        {
                PolygonContainer coincident;
        
                VspBspTraversalData tFrontData;
                VspBspTraversalData tBackData;

                // create new interior node and two leaf nodes
                // or return leaf as it is (if maxCostRatio missed)
                newNode = SubdivideNode(tData, tFrontData, tBackData, coincident);

                if (!newNode->IsLeaf()) //-- continue subdivision
                {
                        // push the children on the stack
                        tStack.push(tFrontData);
                        tStack.push(tBackData);

                        // delete old leaf node
                        DEL_PTR(tData.mNode);   
                }
        }
        
        //-- terminate traversal  
        if (newNode->IsLeaf())
        {
                BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);
        
                // create new view cell for this leaf
                BspViewCell *viewCell = new BspViewCell();
                leaf->SetViewCell(viewCell);
                viewCell->mLeaves.push_back(leaf);

                //-- update pvs
                int conSamp = 0, sampCon = 0;
                AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
                        
                mStat.contributingSamples += conSamp;
                mStat.sampleContributions += sampCon;
                
                EvaluateLeafStats(tData);
        }
        
                
        //-- cleanup
        DEL_PTR(tData.mPolygons);
        DEL_PTR(tData.mRays);
        DEL_PTR(tData.mGeometry);

        return newNode;
}

BspNode *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
                                                                   VspBspTraversalData &frontData,
                                                                   VspBspTraversalData &backData,
                                                                   PolygonContainer &coincident)
{
        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);

        int maxCostMisses = tData.mMaxCostMisses;

        // select subdivision plane
        Plane3 splitPlane;
        if (!SelectPlane(splitPlane, leaf, tData))
        {
                ++ maxCostMisses;

                if (maxCostMisses >= mTermMissTolerance)
                {
                        // terminate branch because of max cost
                        ++ mStat.maxCostNodes; 
            return leaf;
                }
        }

        mStat.nodes += 2;

        //-- subdivide further
        BspInterior *interior = new BspInterior(splitPlane); 

#ifdef _DEBUG
        Debug << interior << endl;
#endif
        
        //-- the front and back traversal data is filled with the new values
        frontData.mPolygons = new PolygonContainer();
        frontData.mDepth = tData.mDepth + 1;
        frontData.mRays = new RayInfoContainer();
        frontData.mGeometry = new BspNodeGeometry();

        backData.mPolygons = new PolygonContainer();
        backData.mDepth = tData.mDepth + 1;
        backData.mRays = new RayInfoContainer();
        backData.mGeometry = new BspNodeGeometry();

        // subdivide rays
        SplitRays(interior->GetPlane(), 
                          *tData.mRays, 
                          *frontData.mRays, 
                          *backData.mRays);
        
        // subdivide polygons
        mStat.splits += SplitPolygons(interior->GetPlane(),
                                                                  *tData.mPolygons, 
                                          *frontData.mPolygons, 
                                                                  *backData.mPolygons,
                                                                  coincident);

        
        // how often was max cost ratio missed in this branch?
        frontData.mMaxCostMisses = maxCostMisses;
        backData.mMaxCostMisses = maxCostMisses;

        // compute pvs
        frontData.mPvs = ComputePvsSize(*frontData.mRays);
        backData.mPvs = ComputePvsSize(*backData.mRays);

        // split geometry and compute area
        if (1)
        {
                tData.mGeometry->SplitGeometry(*frontData.mGeometry, 
                                                                           *backData.mGeometry, 
                                                                           interior->GetPlane(),
                                                                           mBox,
                                                                           mEpsilon);
        
                frontData.mArea = frontData.mGeometry->GetArea();
                backData.mArea = backData.mGeometry->GetArea();
        }

        // compute accumulated ray length
        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);

        //-- create front and back leaf

        BspInterior *parent = leaf->GetParent();

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

        // and setup child links
        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
        
        frontData.mNode = interior->GetFront();
        backData.mNode = interior->GetBack();

        //DEL_PTR(leaf);
        return interior;
}

void VspBspTree::AddToPvs(BspLeaf *leaf,
                                                  const RayInfoContainer &rays, 
                                                  int &sampleContributions,
                                                  int &contributingSamples)
{
        sampleContributions = 0;
        contributingSamples = 0;

    RayInfoContainer::const_iterator it, it_end = rays.end();

        ViewCell *vc = leaf->GetViewCell();

        // add contributions from samples to the PVS
        for (it = rays.begin(); it != it_end; ++ it)
        {
                int contribution = 0;
                VssRay *ray = (*it).mRay;
                        
                if (ray->mTerminationObject)
                        contribution += vc->GetPvs().AddSample(ray->mTerminationObject);
                
                if (ray->mOriginObject)
                        contribution += vc->GetPvs().AddSample(ray->mOriginObject);

                if (contribution)
                {
                        sampleContributions += contribution;
                        ++ contributingSamples;
                }
                //leaf->mVssRays.push_back(ray);
        }
}

void VspBspTree::SortSplitCandidates(const PolygonContainer &polys, 
                                                                         const int axis, 
                                                                         vector<SortableEntry> &splitCandidates) const
{
        splitCandidates.clear();
  
        const int requestedSize = 2 * (int)polys.size();
        
        // creates a sorted split candidates array  
        splitCandidates.reserve(requestedSize);
  
        PolygonContainer::const_iterator it, it_end = polys.end();

        AxisAlignedBox3 box;

        // insert all queries 
        for(it = polys.begin(); it != it_end; ++ it)
        {
                box.Initialize();
                box.Include(*(*it));
          
                splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
                splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
        }
  
        stable_sort(splitCandidates.begin(), splitCandidates.end());
}


float VspBspTree::BestCostRatio(const PolygonContainer &polys,
                                                                const AxisAlignedBox3 &box,
                                                                const int axis,
                                                                float &position,
                                                                int &objectsBack,
                                int &objectsFront) const
{
        vector<SortableEntry> splitCandidates;

        SortSplitCandidates(polys, axis, splitCandidates);
        
        // 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
        
        int objectsLeft = 0, objectsRight = (int)polys.size();
        
        float minBox = box.Min(axis);
        float maxBox = box.Max(axis);
        float boxArea = box.SurfaceArea();
  
        float minBand = minBox + mAxisAlignedSplitBorder * (maxBox - minBox);
        float maxBand = minBox + (1.0f - mAxisAlignedSplitBorder) * (maxBox - minBox);
        
        float minSum = 1e20f;
        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();

        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci) 
        {
                switch ((*ci).type) 
                {
                        case SortableEntry::POLY_MIN:
                                ++ objectsLeft;
                                break;
                        case SortableEntry::POLY_MAX:
                            -- objectsRight;
                                break;
                        default:
                                break;
                }
                
                if ((*ci).value > minBand && (*ci).value < maxBand) 
                {
                        AxisAlignedBox3 lbox = box;
                        AxisAlignedBox3 rbox = box;
                        lbox.SetMax(axis, (*ci).value);
                        rbox.SetMin(axis, (*ci).value);

                        float sum = objectsLeft * lbox.SurfaceArea() + 
                                                objectsRight * rbox.SurfaceArea();
      
                        if (sum < minSum) 
                        {
                                minSum = sum;
                                position = (*ci).value;

                                objectsBack = objectsLeft;
                                objectsFront = objectsRight;
                        }
                }
        }
  
        const float oldCost = (float)polys.size();
        const float newCost = mAxisAlignedCtDivCi + minSum / boxArea;
        const float ratio = newCost / oldCost;


#if 0
  Debug << "====================" << endl;
  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
#endif
  return ratio;
}

bool VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
                                        const PolygonContainer &polys) const
{
        AxisAlignedBox3 box;
        box.Initialize();
        
        // create bounding box of region
        Polygon3::IncludeInBox(polys, box);
        
        int objectsBack = 0, objectsFront = 0;
        int axis = 0;
        float costRatio = MAX_FLOAT;
        Vector3 position;

        //-- area subdivision
        for (int i = 0; i < 3; ++ i) 
        {
                float p = 0;
                const float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
                
                if (r < costRatio)
                {
                        costRatio = r;
                        axis = i;
                        position = p;
                }
        }
        
        if (costRatio >= mTermMaxCostRatio)
                return false;

        Vector3 norm(0,0,0); norm[axis] = 1.0f;
        plane = Plane3(norm, position);

        return true;
}

bool VspBspTree::SelectPlane(Plane3 &plane, 
                                                         BspLeaf *leaf, 
                                                         VspBspTraversalData &data)
{
        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
                ((int)data.mRays->size() > mTermMinRaysForAxisAligned))
        {       // TODO: candidates should be rays
                return SelectAxisAlignedPlane(plane, *data.mPolygons); 
        }

        // simplest strategy: just take next polygon
        if (mSplitPlaneStrategy & RANDOM_POLYGON)
        {
        if (!data.mPolygons->empty())
                {
                        const int randIdx = (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
                        Polygon3 *nextPoly = (*data.mPolygons)[randIdx];

                        plane = nextPoly->GetSupportingPlane();
                        return true;
                }
                else
                {
                        //-- choose plane on midpoint of a ray
                        const int candidateIdx = (int)RandomValue(0, (Real)((int)data.mRays->size() - 1));
                                                                        
                        const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
                        const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();

                        const Vector3 pt = (maxPt + minPt) * 0.5;

                        const Vector3 normal = (*data.mRays)[candidateIdx].mRay->GetDir();
                        
                        plane = Plane3(normal, pt);
                        return true;
                }

                return false;
        }

        // use heuristics to find appropriate plane
        return SelectPlaneHeuristics(plane, leaf, data); 
}

Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
{       
        const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
        
        const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
        const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();

        const Vector3 pt = (maxPt + minPt) * 0.5;
        const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());

        return Plane3(normal, pt);
}

Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
{       
        Vector3 pt[3];
        
        int idx[3];
        int cmaxT = 0;
        int cminT = 0;
        bool chooseMin = false;

        for (int j = 0; j < 3; ++ j)
        {
                idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));
                        
                if (idx[j] >= (int)rays.size())
                {
                        idx[j] -= (int)rays.size();
                                
                        chooseMin = (cminT < 2);
                }
                else
                        chooseMin = (cmaxT < 2);

                RayInfo rayInf = rays[idx[j]];
                pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
        }       

        return Plane3(pt[0], pt[1], pt[2]);
}

Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
{       
        Vector3 pt[3];
        
        int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
        int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));

        // check if rays different
        if (idx1 == idx2)
                idx2 = (idx2 + 1) % (int)rays.size();

        const RayInfo ray1 = rays[idx1];
        const RayInfo ray2 = rays[idx2];

        // normal vector of the plane parallel to both lines
        const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));

        // vector from line 1 to line 2
        const Vector3 vd = (ray2.ExtrapOrigin() - ray1.ExtrapOrigin());
        
        // project vector on normal to get distance
        const float dist = DotProd(vd, norm);

        // point on plane lies halfway between the two planes
        const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;

        return Plane3(norm, planePt);
}

bool VspBspTree::SelectPlaneHeuristics(Plane3 &bestPlane, 
                                                                           BspLeaf *leaf, 
                                                                           VspBspTraversalData &data)
{
        float lowestCost = MAX_FLOAT;
        // intermediate plane
        Plane3 plane;

        const int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
        int maxIdx = (int)data.mPolygons->size();
        
        for (int i = 0; i < limit; ++ i)
        {
                // assure that no index is taken twice
                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
                //Debug << "current Idx: " << maxIdx << " cand idx " << candidateIdx << endl;
                
                Polygon3 *poly = (*data.mPolygons)[candidateIdx];

                // swap candidate to the end to avoid testing same plane
                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
        
                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];

                // evaluate current candidate
                const float candidateCost = 
                        SplitPlaneCost(poly->GetSupportingPlane(), data);

                if (candidateCost < lowestCost)
                {
                        bestPlane = poly->GetSupportingPlane();
                        lowestCost = candidateCost;
                }
        }
        
        //-- choose candidate planes extracted from rays
        //-- different methods are available
        for (int i = 0; i < mMaxRayCandidates; ++ i)
        {
                plane = ChooseCandidatePlane3(*data.mRays);
                const float candidateCost = SplitPlaneCost(plane, data);

                if (candidateCost < lowestCost)
                {
                        bestPlane = plane;      
                        lowestCost = candidateCost;
                }
        }

#ifdef _DEBUG
        Debug << "plane lowest cost: " << lowestCost << endl;
#endif

        // cost ratio miss
        if (lowestCost > mTermMaxCostRatio)
                return false;

        return true;
}


inline void VspBspTree::GenerateUniqueIdsForPvs()
{
        Intersectable::NewMail(); sBackId = ViewCell::sMailId;
        Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
        Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
}

float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
                                                                 const VspBspTraversalData &data)
{
        float cost = 0;

        float sumBalancedRays = 0;
        float sumRaySplits = 0;
        
        int frontPvs = 0;
        int backPvs = 0;

        // probability that view point lies in child
        float pOverall = 0;
        float pFront = 0;
        float pBack = 0;

        const bool pvsUseLen = false;

        if (mSplitPlaneStrategy & PVS)
        {
                // create unique ids for pvs heuristics
                GenerateUniqueIdsForPvs();
        
                if (mPvsUseArea) // use front and back cell areas to approximate volume
                {       
                        // construct child geometry with regard to the candidate split plane
                        BspNodeGeometry frontCell;
                        BspNodeGeometry backCell;
                
                        data.mGeometry->SplitGeometry(frontCell, 
                                                                                  backCell, 
                                                                                  candidatePlane,
                                                                                  mBox,
                                                                                  mEpsilon);
                
                        pFront = frontCell.GetArea();
                        pBack = backCell.GetArea();

                        pOverall = data.mArea;
                }
        }
                
        int limit;
        bool useRand;

        // choose test polyongs randomly if over threshold
        if ((int)data.mRays->size() > mMaxTests)
        {
                useRand = true;
                limit = mMaxTests;
        }
        else
        {
                useRand = false;
                limit = (int)data.mRays->size();
        }

        for (int i = 0; i < limit; ++ i)
        {
                const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;
                RayInfo rayInf = (*data.mRays)[testIdx];

                VssRay *ray = rayInf.mRay;
                float t;
                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);

                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
                {
                        sumBalancedRays += cf;
                }
                
                if (mSplitPlaneStrategy & BALANCED_RAYS)
                {
                        if (cf == 0)
                                ++ sumRaySplits;
                }

                if (mSplitPlaneStrategy & PVS)
                {
                        // in case the ray intersects an object
                        // assure that we only count the object 
                        // once for the front and once for the back side of the plane
                        
                        // add the termination object
                        AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
                        
                        // add the source object
                        AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
                        
                        // use number or length of rays to approximate volume
                        if (!mPvsUseArea)
                        {
                                float len = 1;

                                if (pvsUseLen) // use length of rays
                                        len = rayInf.SqrSegmentLength();
                        
                                pOverall += len;

                                if (cf == 1)
                                        pFront += len;
                                if (cf == -1)
                                        pBack += len;
                                if (cf == 0)
                                {
                                        // use length of rays to approximate volume
                                        if (pvsUseLen) 
                                        {
                                                float newLen = len * 
                                                        (rayInf.GetMaxT() - t) / 
                                                        (rayInf.GetMaxT() - rayInf.GetMinT());
                
                                                if (candidatePlane.Side(rayInf.ExtrapOrigin()) <= 0)
                                                {
                                                        pBack += newLen;
                                                        pFront += len - newLen;
                                                }
                                                else
                                                {
                                                        pFront += newLen;
                                                        pBack += len - newLen;
                                                }
                                        }
                                        else
                                        {
                                                ++ pFront;
                                                ++ pBack;
                                        }
                                }
                        }
                }
        }

        const float raysSize = (float)data.mRays->size() + Limits::Small;
        
        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
                cost += mLeastRaySplitsFactor * sumRaySplits / raysSize;

        if (mSplitPlaneStrategy & BALANCED_RAYS)
                cost += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
        
        // pvs criterium
        if (mSplitPlaneStrategy & PVS)
        {
                const float oldCost = pOverall * (float)data.mPvs + Limits::Small;
                cost += mPvsFactor * (frontPvs * pFront + backPvs * pBack) / oldCost;

                //cost += mPvsFactor * 0.5 * (frontPvs * pFront + backPvs * pBack) / oldCost;
                //Debug << "new cost: " << cost << " over" << frontPvs * pFront + backPvs * pBack << " old cost " << oldCost << endl;
        
                if (0) // give penalty to unbalanced split
                        if (((pFront * 0.2 + Limits::Small) > pBack) || 
                                (pFront < (pBack * 0.2 + Limits::Small)))
                                        cost += 0.5;
        }

#ifdef _DEBUG
        Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
                  << " frontpvs: " << frontPvs << " pFront: " << pFront 
                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
#endif
        
        // normalize cost by sum of linear factors
        return cost / (float)mCostNormalizer;
}

void VspBspTree::AddObjToPvs(Intersectable *obj,
                                                         const int cf,
                                                         int &frontPvs,
                                                         int &backPvs) const
{
        if (!obj)
                return;
        // TODO: does this really belong to no pvs?
        //if (cf == Ray::COINCIDENT) return;

        // object belongs to both PVS
        if (cf >= 0)
        {
                if ((obj->mMailbox != sFrontId) && 
                        (obj->mMailbox != sFrontAndBackId))
                {
                        ++ frontPvs;

                        if (obj->mMailbox == sBackId)
                                obj->mMailbox = sFrontAndBackId;        
                        else
                                obj->mMailbox = sFrontId;                                                               
                }
        }
        
        if (cf <= 0)
        {
                if ((obj->mMailbox != sBackId) &&
                        (obj->mMailbox != sFrontAndBackId))
                {
                        ++ backPvs;

                        if (obj->mMailbox == sFrontId)
                                obj->mMailbox = sFrontAndBackId; 
                        else
                                obj->mMailbox = sBackId;                                
                }
        }
}

void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
{
        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);
 
        while (!nodeStack.empty()) 
        {
                BspNode *node = nodeStack.top();
    
                nodeStack.pop();
    
                if (node->IsLeaf()) 
                {
                        BspLeaf *leaf = (BspLeaf *)node;                
                        leaves.push_back(leaf);
                } 
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

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

AxisAlignedBox3 VspBspTree::GetBoundingBox() const
{
        return mBox;
}

BspNode *VspBspTree::GetRoot() const
{
        return mRoot;
}

void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
{
        // the node became a leaf -> evaluate stats for leafs
        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);

        // store maximal and minimal depth
        if (data.mDepth > mStat.maxDepth)
                mStat.maxDepth = data.mDepth;

        if (data.mDepth < mStat.minDepth)
                mStat.minDepth = data.mDepth;
        
        if (data.mDepth >= mTermMaxDepth)
                ++ mStat.maxDepthNodes;

        if (data.mPvs < mTermMinPvs)
                ++ mStat.minPvsNodes;

        if ((int)data.mRays->size() < mTermMinRays)
                ++ mStat.minRaysNodes;

        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
                ++ mStat.maxRayContribNodes;
        
        if (data.mArea <= mTermMinArea) 
        {
                //Debug << "area: " << data.mArea / mBox.SurfaceArea() << " min area: " << mTermMinArea / mBox.SurfaceArea() << endl;
                ++ mStat.minAreaNodes;
        }
        // accumulate depth to compute average depth
        mStat.accumDepth += data.mDepth;

#ifdef _DEBUG
        Debug << "BSP stats: "
                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
                  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), " 
                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
#endif
}

int VspBspTree::CastRay(Ray &ray)
{
        int hits = 0;
  
        stack<BspRayTraversalData> tStack;
  
        float maxt, mint;

        if (!mBox.GetRaySegment(ray, mint, maxt))
                return 0;

        Intersectable::NewMail();

        Vector3 entp = ray.Extrap(mint);
        Vector3 extp = ray.Extrap(maxt);
  
        BspNode *node = mRoot;
        BspNode *farChild = NULL;
        
        while (1)
        {
                if (!node->IsLeaf()) 
                {
                        BspInterior *in = dynamic_cast<BspInterior *>(node);

                        Plane3 splitPlane = in->GetPlane();
                        const int entSide = splitPlane.Side(entp);
                        const int extSide = splitPlane.Side(extp);

                        if (entSide < 0)
                        {
                                node = in->GetBack();

                                if(extSide <= 0) // plane does not split ray => no far child
                                        continue;
                                        
                                farChild = in->GetFront(); // plane splits ray

                        } else if (entSide > 0)
                        {
                                node = in->GetFront();

                                if (extSide >= 0) // plane does not split ray => no far child
                                        continue;

                                farChild = in->GetBack(); // plane splits ray                   
                        }
                        else // ray and plane are coincident
                        {
                                // WHAT TO DO IN THIS CASE ?
                                //break;
                                node = in->GetFront();
                                continue;
                        }

                        // push data for far child
                        tStack.push(BspRayTraversalData(farChild, extp, maxt));

                        // find intersection of ray segment with plane
                        float t;
                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
                        maxt *= t;
                        
                } else // reached leaf => intersection with view cell
                {
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
      
                        if (!leaf->GetViewCell()->Mailed())
                        {
                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
                                leaf->GetViewCell()->Mail();
                                ++ hits;
                        }
                        
                        //-- fetch the next far child from the stack
                        if (tStack.empty())
                                break;
      
                        entp = extp;
                        mint = maxt; // NOTE: need this?

                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
                                break;

                        BspRayTraversalData &s = tStack.top();

                        node = s.mNode;
                        extp = s.mExitPoint;
                        maxt = s.mMaxT;

                        tStack.pop();
                }
        }

        return hits;
}

bool VspBspTree::Export(const string filename)
{
        Exporter *exporter = Exporter::GetExporter(filename);

        if (exporter) 
        {
                //exporter->ExportVspBspTree(*this);
                return true;
        }       

        return false;
}

void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
{
        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);

        ViewCell::NewMail();

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

                if (node->IsLeaf()) 
                {
                        ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();

                        if (!viewCell->Mailed()) 
                        {
                                viewCell->Mail();
                                viewCells.push_back(viewCell);
                        }
                }
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

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

void VspBspTree::EvaluateViewCellsStats(ViewCellsStatistics &stat) const
{
        stat.Reset();

        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);

        ViewCell::NewMail();

        // exclude root cell
        mRootCell->Mail();

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

                if (node->IsLeaf()) 
                {
                        ++ stat.leaves;

                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();

                        if (!viewCell->Mailed()) 
                        {
                                viewCell->Mail();
                                
                                ++ stat.viewCells;
                                const int pvsSize = viewCell->GetPvs().GetSize();

                stat.pvs += pvsSize;

                                if (pvsSize < 1)
                                        ++ stat.emptyPvs;

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

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

                                if ((int)viewCell->mLeaves.size() > stat.maxLeaves)
                                        stat.maxLeaves = (int)viewCell->mLeaves.size();         
                        }
                }
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

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

BspTreeStatistics &VspBspTree::GetStat()
{
        return mStat;
}

float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
{
        float len = 0;

        RayInfoContainer::const_iterator it, it_end = rays.end();

        for (it = rays.begin(); it != it_end; ++ it)
                len += (*it).SegmentLength();

        return len;
}

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

        while (!rays.empty())
        {
                RayInfo bRay = rays.back();
                rays.pop_back();

                VssRay *ray = bRay.mRay;
                float t;

                        // get classification and receive new t
                const int cf = bRay.ComputeRayIntersection(plane, t);
        
                switch (cf)
                {
                case -1:
                        backRays.push_back(bRay);
                        break;
                case 1:
                        frontRays.push_back(bRay);
                        break;
                case 0: 
                        //-- split ray
                        //--  look if start point behind or in front of plane
                        if (plane.Side(bRay.ExtrapOrigin()) <= 0)
                        {       
                                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()));
                        }
                        break;
                default:
                        Debug << "Should not come here 4" << endl;
                        break;
                }
        }

        return splits;
}

void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
{
        BspNode *lastNode;

        do
        {
                lastNode = n;

                // want to get planes defining geometry of this node => don't take
                // split plane of node itself
                n = n->GetParent();
                
                if (n)
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();

            if (interior->GetFront() != lastNode)
                                halfSpace.ReverseOrientation();

                        halfSpaces.push_back(halfSpace);
                }
        }
        while (n);
}

void VspBspTree::ConstructGeometry(BspNode *n, 
                                                                   BspNodeGeometry &cell) const
{
        PolygonContainer polys;
        ConstructGeometry(n, polys);
        cell.mPolys = polys;
}

void VspBspTree::ConstructGeometry(BspNode *n, 
                                                                   PolygonContainer &cell) const
{
        vector<Plane3> halfSpaces;
        ExtractHalfSpaces(n, halfSpaces);

        PolygonContainer candidatePolys;

        // bounded planes are added to the polygons (reverse polygons
        // as they have to be outfacing
        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
        {
                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
                
                if (p->Valid(mEpsilon))
                {
                        candidatePolys.push_back(p->CreateReversePolygon());
                        DEL_PTR(p);
                }
        }

        // add faces of bounding box (also could be faces of the cell)
        for (int i = 0; i < 6; ++ i)
        {
                VertexContainer vertices;
        
                for (int j = 0; j < 4; ++ j)
                        vertices.push_back(mBox.GetFace(i).mVertices[j]);

                candidatePolys.push_back(new Polygon3(vertices));
        }

        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
        {
                // polygon is split by all other planes
                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
                {
                        if (i == j) // polygon and plane are coincident
                                continue;

                        VertexContainer splitPts;
                        Polygon3 *frontPoly, *backPoly;

                        const int cf = 
                                candidatePolys[i]->ClassifyPlane(halfSpaces[j],
                                                                                                 mEpsilon);
                        
                        switch (cf)
                        {
                                case Polygon3::SPLIT:
                                        frontPoly = new Polygon3();
                                        backPoly = new Polygon3();

                                        candidatePolys[i]->Split(halfSpaces[j], 
                                                                                         *frontPoly, 
                                                                                         *backPoly,
                                                                                         mEpsilon);

                                        DEL_PTR(candidatePolys[i]);

                                        if (frontPoly->Valid(mEpsilon))
                                                candidatePolys[i] = frontPoly;
                                        else
                                                DEL_PTR(frontPoly);

                                        DEL_PTR(backPoly);
                                        break;
                                case Polygon3::BACK_SIDE:
                                        DEL_PTR(candidatePolys[i]);
                                        break;
                                // just take polygon as it is
                                case Polygon3::FRONT_SIDE:
                                case Polygon3::COINCIDENT:
                                default:
                                        break;
                        }
                }
                
                if (candidatePolys[i])
                        cell.push_back(candidatePolys[i]);
        }
}

void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &vcGeom) const
{
        vector<BspLeaf *> leaves = vc->mLeaves;

        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();

        for (it = leaves.begin(); it != it_end; ++ it)
                ConstructGeometry(*it, vcGeom);
}

int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors, 
                                                   const bool onlyUnmailed) const
{
        PolygonContainer cell;

        ConstructGeometry(n, cell);

        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);
                
        // planes needed to verify that we found neighbor leaf.
        vector<Plane3> halfSpaces;
        ExtractHalfSpaces(n, halfSpaces);

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

                if (node->IsLeaf())
                {
            if (node != n && (!onlyUnmailed || !node->Mailed()))
                        {
                                // test all planes of current node if candidate really
                                // is neighbour
                                PolygonContainer neighborCandidate;
                                ConstructGeometry(node, neighborCandidate);
                                
                                bool isAdjacent = true;
                                for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
                                {
                                        const int cf = 
                                                Polygon3::ClassifyPlane(neighborCandidate, 
                                                                                                halfSpaces[i],
                                                                                                mEpsilon);

                                        if (cf == Polygon3::BACK_SIDE)
                                                isAdjacent = false;
                                }

                                if (isAdjacent)
                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));

                                CLEAR_CONTAINER(neighborCandidate);
                        }
                }
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
        
                        const int cf = Polygon3::ClassifyPlane(cell, 
                                                                                                   interior->GetPlane(), 
                                                                                                   mEpsilon);

                        if (cf == Polygon3::FRONT_SIDE)
                                nodeStack.push(interior->GetFront());
                        else
                                if (cf == Polygon3::BACK_SIDE)
                                        nodeStack.push(interior->GetBack());
                                else 
                                {
                                        // random decision
                                        nodeStack.push(interior->GetBack());
                                        nodeStack.push(interior->GetFront());
                                }
                }
        }
        
        CLEAR_CONTAINER(cell);
        return (int)neighbors.size();
}

BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
{
    stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);
        
        int mask = rand();
  
        while (!nodeStack.empty()) 
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
          
                if (node->IsLeaf()) 
                {
                        return dynamic_cast<BspLeaf *>(node);
                } 
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
                        
                        BspNode *next;
        
                        PolygonContainer cell;

                        // todo: not very efficient: constructs full cell everytime
                        ConstructGeometry(interior, cell);

                        const int cf = Polygon3::ClassifyPlane(cell, halfspace, mEpsilon);

                        if (cf == Polygon3::BACK_SIDE)
                                next = interior->GetFront();
                        else
                                if (cf == Polygon3::FRONT_SIDE)
                                        next = interior->GetFront();
                        else 
                        {
                                // random decision
                                if (mask & 1)
                                        next = interior->GetBack();
                                else
                                        next = interior->GetFront();
                                mask = mask >> 1;
                        }

                        nodeStack.push(next);
                }
        }
        
        return NULL;
}

BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
{
        stack<BspNode *> nodeStack;
        
        nodeStack.push(mRoot);

        int mask = rand();
        
        while (!nodeStack.empty()) 
        {
                BspNode *node = nodeStack.top();
                nodeStack.pop();
                
                if (node->IsLeaf()) 
                {
                        if ( (!onlyUnmailed || !node->Mailed()) )
                                return dynamic_cast<BspLeaf *>(node);
                }
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

                        // random decision
                        if (mask & 1)
                                nodeStack.push(interior->GetBack());
                        else
                                nodeStack.push(interior->GetFront());

                        mask = mask >> 1;
                }
        }
        
        return NULL;
}

int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
{
        int pvsSize = 0;

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

        Intersectable::NewMail();

        for (rit = rays.begin(); rit != rays.end(); ++ rit)
        {
                VssRay *ray = (*rit).mRay;
                
                if (ray->mOriginObject)
                {
                        if (!ray->mOriginObject->Mailed())
                        {
                                ray->mOriginObject->Mail();
                                ++ pvsSize;
                        }
                }
                if (ray->mTerminationObject)
                {
                        if (!ray->mTerminationObject->Mailed())
                        {
                                ray->mTerminationObject->Mail();
                                ++ pvsSize;
                        }
                }
        }

        return pvsSize;
}

float VspBspTree::GetEpsilon() const
{
        return mEpsilon;
}

BspViewCell *VspBspTree::GetRootCell() const
{
        return mRootCell;
}

int VspBspTree::SplitPolygons(const Plane3 &plane,
                                                          PolygonContainer &polys, 
                                                          PolygonContainer &frontPolys, 
                                                          PolygonContainer &backPolys, 
                                                          PolygonContainer &coincident) const
{
        int splits = 0;

        while (!polys.empty())
        {
                Polygon3 *poly = polys.back();
                polys.pop_back();

                // classify polygon
                const int cf = poly->ClassifyPlane(plane, mEpsilon);

                switch (cf)
                {
                        case Polygon3::COINCIDENT:
                                coincident.push_back(poly);
                                break;                  
                        case Polygon3::FRONT_SIDE:      
                                frontPolys.push_back(poly);
                                break;
                        case Polygon3::BACK_SIDE:
                                backPolys.push_back(poly);
                                break;
                        case Polygon3::SPLIT:
                                backPolys.push_back(poly);
                                frontPolys.push_back(poly);
                                ++ splits;
                                break;
                        default:
                Debug << "SHOULD NEVER COME HERE\n";
                                break;
                }
        }

        return splits;
}


int VspBspTree::CastLineSegment(const Vector3 &origin,
                                                                const Vector3 &termination,
                                                                vector<ViewCell *> &viewcells)
{
        int hits = 0;
        stack<BspRayTraversalData> tStack;
  
        float mint = 0.0f, maxt = 1.0f;
  
        Intersectable::NewMail();
  
        Vector3 entp = origin;
        Vector3 extp = termination;
  
        BspNode *node = mRoot;
        BspNode *farChild = NULL;
  
        while (1) 
        {
                if (!node->IsLeaf())  
                {
                        BspInterior *in = dynamic_cast<BspInterior *>(node);
          
                        Plane3 splitPlane = in->GetPlane();
                        const int entSide = splitPlane.Side(entp);
                        const int extSide = splitPlane.Side(extp);
          
                        if (entSide < 0) 
                        {
                                node = in->GetBack();
                
                                if(extSide <= 0) // plane does not split ray => no far child
                                        continue;
                
                                farChild = in->GetFront(); // plane splits ray
                        } else if (entSide > 0) 
                        {
                                node = in->GetFront();
                  
                                if (extSide >= 0) // plane does not split ray => no far child
                                        continue;
                  
                                farChild = in->GetBack(); // plane splits ray                   
                        }
                        else // ray and plane are coincident
                        {
                                // WHAT TO DO IN THIS CASE ?
                                //break;
                                node = in->GetFront();
                                continue;
                        }
          
                        // push data for far child
                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
          
                        // find intersection of ray segment with plane
                        float t;
                        extp = splitPlane.FindIntersection(origin, extp, &t);
                        maxt *= t;
          
                } else 
                {
                        // reached leaf => intersection with view cell
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
          
                        if (!leaf->GetViewCell()->Mailed()) 
                        {
                                viewcells.push_back(leaf->GetViewCell());
                                leaf->GetViewCell()->Mail();
                                ++ hits;
                        }
          
                        //-- fetch the next far child from the stack
                        if (tStack.empty())
                                break;
      
                        entp = extp;
                        mint = maxt; // NOTE: need this?
          
                        
                        BspRayTraversalData &s = tStack.top();
            
                        node = s.mNode;
                        extp = s.mExitPoint;
                        maxt = s.mMaxT;
          
                        tStack.pop();
                }
        }
        return hits;
}