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

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

#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 "Exporter.h"
#include "Plane3.h"
#include "ViewCellBsp.h"
#include "ViewCellsManager.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;

float BspMergeCandidate::sOverallCost = Limits::Small;

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

VspBspTree::VspBspTree():
mRoot(NULL),
mPvsUseArea(true),
mCostNormalizer(Limits::Small),
mViewCellsManager(NULL),
mStoreRays(true),
mOnlyDrivingAxis(false)
{
        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);

        //-- max cost ratio for early tree termination
        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);

        //-- partition criteria
        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);

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

        // maximum and minimum number of view cells
        environment->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
        environment->GetIntValue("VspBspTree.PostProcess.minViewCells", mMergeMinViewCells);
        environment->GetFloatValue("VspBspTree.PostProcess.maxCostRatio", mMergeMaxCostRatio);


        //-- debug output
        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";
        }

        mSplitCandidates = new vector<SortableEntry>;

        Debug << endl;
}


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


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

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
{
        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))
        {
                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 and create new view cell
        if (newNode->IsLeaf())
        {
                BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);

                // create new view cell for this leaf
                BspViewCell *viewCell = new BspViewCell();
                leaf->SetViewCell(viewCell);

                if (mStoreRays)
                {
                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
                        for (it = tData.mRays->begin(); it != it_end; ++ it)
                                leaf->mVssRays.push_back((*it).mRay);
                }

                viewCell->mLeaves.push_back(leaf);
                viewCell->SetArea(tData.mArea);

                //-- update pvs
                int conSamp = 0, sampCon = 0;
                AddToPvs(leaf, *tData.mRays, conSamp, sampCon);

                mStat.contributingSamples += conSamp;
                mStat.sampleContributions += sampCon;

                EvaluateLeafStats(tData);
        }


        //-- cleanup
        tData.Clear();

        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 RayInfoContainer &rays, const int axis)
{
        mSplitCandidates->clear();

        int requestedSize = 2 * (int)(rays.size());
        // creates a sorted split candidates array
        if (mSplitCandidates->capacity() > 500000 &&
                requestedSize < (int)(mSplitCandidates->capacity()/10) )
        {
        delete mSplitCandidates;
                mSplitCandidates = new vector<SortableEntry>;
        }

        mSplitCandidates->reserve(requestedSize);

        // insert all queries
        for(RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
        {
                bool positive = (*ri).mRay->HasPosDir(axis);
                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
                                                                                                  (*ri).ExtrapOrigin(axis), (*ri).mRay));

                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
                                                                                                  (*ri).ExtrapTermination(axis), (*ri).mRay));
        }

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

float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
                                                                                  const AxisAlignedBox3 &box,
                                                                                  const int pvsSize,
                                                                                  const int &axis,
                                          float &position)
{
        int raysBack;
        int raysFront;
        int pvsBack;
        int pvsFront;

        SortSplitCandidates(rays, axis);

        // go through the lists, count the number of objects left and right
        // and evaluate the following cost funcion:
        // C = ct_div_ci  + (ql*rl + qr*rr)/queries

        int rl = 0, rr = (int)rays.size();
        int pl = 0, pr = pvsSize;

        float minBox = box.Min(axis);
        float maxBox = box.Max(axis);
        float sizeBox = maxBox - minBox;

        float minBand = minBox + 0.1f*(maxBox - minBox);
        float maxBand = minBox + 0.9f*(maxBox - minBox);

        float sum = rr*sizeBox;
        float minSum = 1e20f;

        Intersectable::NewMail();

        // set all object as belonging to the fron pvs
        for(RayInfoContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
        {
                if ((*ri).mRay->IsActive())
                {
                        Intersectable *object = (*ri).mRay->mTerminationObject;

                        if (object)
                        {
                                if (!object->Mailed())
                                {
                                        object->Mail();
                                        object->mCounter = 1;
                                }
                                else
                                        ++ object->mCounter;
                        }
                }
        }

        Intersectable::NewMail();

        for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
                ci < mSplitCandidates->end(); ++ ci)
        {
                VssRay *ray;

                switch ((*ci).type)
                {
                        case SortableEntry::ERayMin:
                                {
                                        ++ rl;
                                        ray = (VssRay *) (*ci).ray;

                                        Intersectable *object = ray->mTerminationObject;

                                        if (object && !object->Mailed())
                                        {
                                                object->Mail();
                                                ++ pl;
                                        }
                                        break;
                                }
                        case SortableEntry::ERayMax:
                                {
                                        -- rr;
                                        ray = (VssRay *) (*ci).ray;

                                        Intersectable *object = ray->mTerminationObject;

                                        if (object)
                                        {
                                                if (-- object->mCounter == 0)
                                                        -- pr;
                                        }

                                        break;
                                }
                }

                // Note: sufficient to compare size of bounding boxes of front and back side?
                if ((*ci).value > minBand && (*ci).value < maxBand)
                {
                        sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);

                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
                        // cout<<"cost= "<<sum<<endl;

                        if (sum < minSum)
                        {
                                minSum = sum;
                                position = (*ci).value;

                                raysBack = rl;
                                raysFront = rr;

                                pvsBack = pl;
                                pvsFront = pr;

                        }
                }
        }

        float oldCost = (float)pvsSize;
        float newCost = mCtDivCi + minSum / sizeBox;
        float ratio = newCost / oldCost;

        //Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
        //     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;

        return ratio;
}


float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
                                                                                 const VspBspTraversalData &tData)
{
        AxisAlignedBox3 box;
        box.Initialize();

        // create bounding box of region
        RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();

        for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
                box.Include((*ri).ExtrapTermination());

        const bool useCostHeuristics = false;

        //-- regular split
        float nPosition[3];
        float nCostRatio[3];
        int bestAxis = -1;

        const int sAxis = box.Size().DrivingAxis();

        for (int axis = 0; axis < 3; ++ axis)
        {
                if (!mOnlyDrivingAxis || axis == sAxis)
                {
                        if (!useCostHeuristics)
                        {
                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
                                Vector3 normal(0,0,0); normal[axis] = 1;

                                nCostRatio[axis] = SplitPlaneCost(Plane3(normal, nPosition[axis]), tData);
                        }
                        else
                        {
                                nCostRatio[axis] =
                                        BestCostRatioHeuristics(*tData.mRays,
                                                                                    box,
                                                                                        tData.mPvs,
                                                                                        axis,
                                                                                        nPosition[axis]);
                        }

                        if (bestAxis == -1)
                        {
                                bestAxis = axis;
                        }

                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
                        {
                                /*Debug << "pvs front " << nPvsBack[axis]
                                          << " pvs back " << nPvsFront[axis]
                                          << " overall pvs " << leaf->GetPvsSize() << endl;*/
                                bestAxis = axis;
                        }

                }
        }

        //-- assign best axis
        Vector3 normal(0,0,0); normal[bestAxis] = 1;
        plane = Plane3(normal, nPosition[bestAxis]);

        return nCostRatio[bestAxis];
}


bool VspBspTree::SelectPlane(Plane3 &plane,
                                                         BspLeaf *leaf,
                                                         VspBspTraversalData &data)
{
        // 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();

        float candidateCost;

        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
                candidateCost = SplitPlaneCost(poly->GetSupportingPlane(), data);

                if (candidateCost < lowestCost)
                {
                        bestPlane = poly->GetSupportingPlane();
                        lowestCost = candidateCost;
                }
        }

#if 0
        //-- choose candidate planes extracted from rays
        //-- different methods are available
        for (int i = 0; i < mMaxRayCandidates; ++ i)
        {
                plane = ChooseCandidatePlane3(*data.mRays);
                candidateCost = SplitPlaneCost(plane, data);

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

        // axis aligned splits
        candidateCost = SelectAxisAlignedPlane(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) const
{
        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());
                }
        }
}


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

int VspBspTree::TreeDistance(BspNode *n1, BspNode *n2)
{
        std::deque<BspNode *> path1;
        BspNode *p1 = n1;

        // create path from node 1 to root
        while (p1)
        {
                if (p1 == n2) // second node on path
                        return (int)path1.size();

                path1.push_front(p1);
                p1 = p1->GetParent();
        }

        int depth = n2->GetDepth();
        int d = depth;

        BspNode *p2 = n2;

        // compare with same depth
        while (1)
        {
                if ((d < (int)path1.size()) && (p2 == path1[d]))
                        return (depth - d) + ((int)path1.size() - 1 - d);

                -- d;
                p2 = p2->GetParent();
        }

        return 0; // never come here
}

BspNode *VspBspTree::CollapseTree(BspNode *node)
{
    if (node->IsLeaf())
                return node;

        BspInterior *interior = dynamic_cast<BspInterior *>(node);

        BspNode *front = CollapseTree(interior->GetFront());
        BspNode *back = CollapseTree(interior->GetBack());

        if (front->IsLeaf() && back->IsLeaf())
        {
                BspLeaf *frontLeaf = dynamic_cast<BspLeaf *>(front);
                BspLeaf *backLeaf = dynamic_cast<BspLeaf *>(back);

                //-- collapse tree
                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
                {
                        BspViewCell *vc = frontLeaf->GetViewCell();

                        BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);

                        // replace a link from node's parent
                        if (leaf->GetParent())
                                leaf->GetParent()->ReplaceChildLink(node, leaf);

                        delete interior;

                        return leaf;
                }
        }

        return node;
}


void VspBspTree::RepairVcLeafLists()
{
        // list not valid anymore => clear
        stack<BspNode *> nodeStack;
        nodeStack.push(mRoot);

        ViewCell::NewMail();

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

                if (node->IsLeaf())
                {
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);

                        BspViewCell *viewCell = leaf->GetViewCell();

                        if (!viewCell->Mailed())
                        {
                                viewCell->mLeaves.clear();
                                viewCell->Mail();
                        }

                        viewCell->mLeaves.push_back(leaf);
                }
                else
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

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


int VspBspTree::MergeLeaves()
{
        vector<BspLeaf *> leaves;
        priority_queue<BspMergeCandidate> mergeQueue;

        // collect the leaves, e.g., the "voxels" that will build the view cells
        CollectLeaves(leaves);

        int vcSize = (int)leaves.size();
        int savedVcSize = vcSize;

        BspLeaf::NewMail();

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

        // find merge candidates and push them into queue
        for (it = leaves.begin(); it != it_end; ++ it)
        {
                BspLeaf *leaf = *it;

                // no leaf is part of two merge candidates
                if (!leaf->Mailed())
                {
                        leaf->Mail();

                        vector<BspLeaf *> neighbors;
                        FindNeighbors(leaf, neighbors, true);

                        vector<BspLeaf *>::const_iterator nit,
                                nit_end = neighbors.end();

                        for (nit = neighbors.begin(); nit != nit_end; ++ nit)
                        {
                                BspMergeCandidate mc = BspMergeCandidate(leaf, *nit);
                                mergeQueue.push(mc);

                                BspMergeCandidate::sOverallCost += mc.GetLeaf1Cost();
                                BspMergeCandidate::sOverallCost += mc.GetLeaf2Cost();
                        }
                }
        }

        int merged = 0;
        int mergedSiblings = 0;
        int accTreeDist = 0;
        int maxTreeDist = 0;
        const bool mergeStats = true;


        //-- use priority queue to merge leaf pairs
        while (!mergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
                   (mergeQueue.top().GetMergeCost() <
                    mMergeMaxCostRatio * BspMergeCandidate::sOverallCost))
        {
                //Debug << "mergecost: " << mergeQueue.top().GetMergeCost() / BspMergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
                BspMergeCandidate mc = mergeQueue.top();
                mergeQueue.pop();

                // both view cells equal!
                if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
                        continue;

                if (mc.Valid())
                {
                        ViewCell::NewMail();
                        MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
                        -- vcSize;
                        // increase absolute merge cost
                        BspMergeCandidate::sOverallCost += mc.GetMergeCost();

                        //-- stats
                        ++ merged;

                        if (mc.GetLeaf1()->IsSibling(mc.GetLeaf2()))
                                ++ mergedSiblings;

                        if (mergeStats)
                        {
                                const int dist = TreeDistance(mc.GetLeaf1(), mc.GetLeaf2());

                                if (dist > maxTreeDist)
                                        maxTreeDist = dist;

                                accTreeDist += dist;
                        }
                }
                // merge candidate not valid, because one of the leaves was already
                // merged with another one
                else
                {
                        // validate and reinsert into queue
                        mc.SetValid();
                        mergeQueue.push(mc);
                        //Debug << "validate " << mc.GetMergeCost() << endl;
                }
        }

        // collapse tree according to view cell partition
        CollapseTree(mRoot);
        // revalidate leaves
        RepairVcLeafLists();

        Debug << "Merged " << merged << " nodes of " << savedVcSize
                  << " (merged " << mergedSiblings << " siblings)" << endl;

        if (mergeStats)
        {
                Debug << "maximal tree distance: " << maxTreeDist << endl;
                Debug << "Avg tree distance: " << (float)accTreeDist / (float)merged << endl;
        }

        //TODO: should return sample contributions
        return merged;
}


bool VspBspTree::MergeViewCells(BspLeaf *l1, BspLeaf *l2)
{
        //-- change pointer to view cells of all leaves associated
        //-- with the previous view cells
        BspViewCell *fVc = l1->GetViewCell();
        BspViewCell *bVc = l2->GetViewCell();

        BspViewCell *vc = dynamic_cast<BspViewCell *>(
                mViewCellsManager->MergeViewCells(*fVc, *bVc));

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

        // set new size of view cell
        vc->SetArea(fVc->GetArea() + bVc->GetArea());

        vector<BspLeaf *> fLeaves = fVc->mLeaves;
        vector<BspLeaf *> bLeaves = bVc->mLeaves;

        vector<BspLeaf *>::const_iterator it;

        //-- change view cells of all the other leaves the view cell belongs to
        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
        {
                (*it)->SetViewCell(vc);
                vc->mLeaves.push_back(*it);
        }

        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
        {
                (*it)->SetViewCell(vc);
                vc->mLeaves.push_back(*it);
        }

        vc->Mail();

        // clean up old view cells
        DEL_PTR(fVc);
        DEL_PTR(bVc);

        return true;
}


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

/************************************************************************/
/*                BspMergeCandidate implementation                      */
/************************************************************************/


BspMergeCandidate::BspMergeCandidate(BspLeaf *l1, BspLeaf *l2):
mMergeCost(0),
mLeaf1(l1),
mLeaf2(l2),
mLeaf1Id(l1->GetViewCell()->mMailbox),
mLeaf2Id(l2->GetViewCell()->mMailbox)
{
        EvalMergeCost();
}

float BspMergeCandidate::GetLeaf1Cost() const
{
        BspViewCell *vc = mLeaf1->GetViewCell();
        return vc->GetPvs().GetSize() * vc->GetArea();
}

float BspMergeCandidate::GetLeaf2Cost() const
{
        BspViewCell *vc = mLeaf2->GetViewCell();
        return vc->GetPvs().GetSize() * vc->GetVolume();
}

void BspMergeCandidate::EvalMergeCost()
{
        //-- compute pvs difference
        BspViewCell *vc1 = mLeaf1->GetViewCell();
        BspViewCell *vc2 = mLeaf2->GetViewCell();

        const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
        const int vcPvs = diff1 + vc1->GetPvs().GetSize();

        //-- compute ratio of old cost
        //-- (added size of left and right view cell times pvs size)
        //-- to new rendering cost (size of merged view cell times pvs size)
        const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();

        const float newCost =
                (float)vcPvs * (vc1->GetArea() + vc2->GetArea());

        mMergeCost = newCost - oldCost;

//      if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
//              mMergeCost += 1.0;
}

void BspMergeCandidate::SetLeaf1(BspLeaf *l)
{
        mLeaf1 = l;
}

void BspMergeCandidate::SetLeaf2(BspLeaf *l)
{
        mLeaf2 = l;
}

BspLeaf *BspMergeCandidate::GetLeaf1()
{
        return mLeaf1;
}

BspLeaf *BspMergeCandidate::GetLeaf2()
{
        return mLeaf2;
}

bool BspMergeCandidate::Valid() const
{
        return
                (mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
                (mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
}

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

void BspMergeCandidate::SetValid()
{
        mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
        mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;

        EvalMergeCost();
}