source: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp @ 367

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

int BspTree::sTermMaxPolygons = 10;
int BspTree::sTermMaxDepth = 20;
int BspTree::sMaxCandidates = 10;
int BspTree::sSplitPlaneStrategy = BALANCED_POLYS; 
int BspTree::sConstructionMethod = FROM_INPUT_VIEW_CELLS;
int BspTree::sTermMaxPolysForAxisAligned = 50;
int BspTree::sTermMaxObjectsForAxisAligned = 50;
int BspTree::sTermMaxRaysForAxisAligned = -1;
int BspTree::sTermMaxRays = -1;
int BspTree::sMinPvsDif = 10;

float BspTree::sCt_div_ci = 1.0f;
float BspTree::sSplitBorder = 1.0f;
float BspTree::sMaxCostRatio = 0.9f;

//-- factors for bsp tree split plane heuristics

float BspTree::sLeastSplitsFactor = 1.0f;
float BspTree::sBalancedPolysFactor = 1.0f;
float BspTree::sBalancedViewCellsFactor = 1.0f;

// NOTE:  very important criterium for 2.5d scenes
float BspTree::sVerticalSplitsFactor = 1.0f;
float BspTree::sLargestPolyAreaFactor = 1.0f;
float BspTree::sBlockedRaysFactor = 1.0f;
float BspTree::sLeastRaySplitsFactor = 1.0f;
float BspTree::sBalancedRaysFactor = 1.0f;
float BspTree::sPvsFactor = 1.0f;

bool BspTree::sStoreSplitPolys = false;

int BspNode::mailID = 1;

/** Evaluates split plane classification with respect to the plane's 
        contribution for a minimum number splits in the tree.
*/
float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
/** Evaluates split plane classification with respect to the plane's 
        contribution for a balanced tree.
*/
float BspTree::sBalancedPolysTable[] = {1, -1, 0, 0};

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

/****************************************************************/
/*                  class BspNode implementation                */
/****************************************************************/

BspNode::BspNode(): 
mParent(NULL), mPolygons(NULL)
{}

BspNode::BspNode(BspInterior *parent): 
mParent(parent), mPolygons(NULL)
{}

BspNode::~BspNode()
{
        if (mPolygons)
        {
                CLEAR_CONTAINER(*mPolygons);
                DEL_PTR(mPolygons);
        }
}

bool BspNode::IsRoot() const 
{
        return mParent == NULL;
}

BspInterior *BspNode::GetParent() 
{ 
        return mParent; 
}

void BspNode::SetParent(BspInterior *parent)
{
        mParent = parent;
}

PolygonContainer *BspNode::GetPolygons()
{
        if (!mPolygons)
                mPolygons = new PolygonContainer();

        return mPolygons;
}

void BspNode::ProcessPolygons(PolygonContainer *polys, bool storePolys)
{
        if (storePolys)
        {
                while (!polys->empty())
                {
                        GetPolygons()->push_back(polys->back());
                        polys->pop_back();
                }
        }
        else CLEAR_CONTAINER(*polys);
}


/****************************************************************/
/*              class BspInterior implementation                */
/****************************************************************/


BspInterior::BspInterior(const Plane3 &plane): 
mPlane(plane), mFront(NULL), mBack(NULL)
{}

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

BspNode *BspInterior::GetBack() 
{
        return mBack;
}

BspNode *BspInterior::GetFront() 
{
        return mFront;
}

Plane3 *BspInterior::GetPlane()
{
        return &mPlane;
}

void BspInterior::ReplaceChildLink(BspNode *oldChild, BspNode *newChild) 
{
        if (mBack == oldChild)
                mBack = newChild;
        else
                mFront = newChild;
}

void BspInterior::SetupChildLinks(BspNode *b, BspNode *f) 
{
    mBack = b;
    mFront = f;
}

void BspInterior::ProcessPolygon(Polygon3 **poly, const bool storePolys)
{
        if (storePolys) 
                GetPolygons()->push_back(*poly);
        else
                DEL_PTR(*poly);
}

int BspInterior::SplitPolygons(PolygonContainer &polys, 
                                                           PolygonContainer &frontPolys, 
                                                           PolygonContainer &backPolys, 
                                                           PolygonContainer &coincident,
                                                           const bool storePolys)
{
        Polygon3 *splitPoly = NULL;

        int splits = 0;

#ifdef _Debug
        Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
#endif
        while (!polys.empty())
        {
                Polygon3 *poly = polys.back();
                polys.pop_back();

                //Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";

                // classify polygon
                const int cf = poly->ClassifyPlane(mPlane);

                Polygon3 *front_piece = NULL;
                Polygon3 *back_piece = NULL;
        
                VertexContainer splitVertices;

                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:
                                front_piece = new Polygon3(poly->mParent);
                                back_piece = new Polygon3(poly->mParent);

                                //-- split polygon into front and back part
                                poly->Split(mPlane, 
                                                        *front_piece, 
                                                        *back_piece, 
                                                        splitVertices);
                                        
                                ++ splits; // increase number of splits

                                //-- inherit rays from parent polygon for blocked ray criterium
                                poly->InheritRays(*front_piece, *back_piece);
                                //Debug << "p: " << poly->mPiercingRays.size() << " f: " << front_piece->mPiercingRays.size() << " b: " << back_piece->mPiercingRays.size() << endl;
                                
                                // check if polygons still valid 
                                if (front_piece->Valid())
                                        frontPolys.push_back(front_piece);
                                else
                                        DEL_PTR(front_piece);
                                
                                if (back_piece->Valid())
                                        backPolys.push_back(back_piece);
                                else                            
                                        DEL_PTR(back_piece);
                                
#ifdef _DEBUG
                                Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
#endif
                                ProcessPolygon(&poly, storePolys);                      
                                break;
                        default:
                Debug << "SHOULD NEVER COME HERE\n";
                                break;
                }
        }

        return splits;
}

/****************************************************************/
/*                  class BspLeaf implementation                */
/****************************************************************/
BspLeaf::BspLeaf(): mViewCell(NULL)
{
}

BspLeaf::BspLeaf(BspViewCell *viewCell): 
mViewCell(viewCell)
{
}

BspLeaf::BspLeaf(BspInterior *parent): 
BspNode(parent), mViewCell(NULL)
{}

BspLeaf::BspLeaf(BspInterior *parent, BspViewCell *viewCell): 
BspNode(parent), mViewCell(viewCell)
{
}

BspViewCell *BspLeaf::GetViewCell() const
{
        return mViewCell;
}

void BspLeaf::SetViewCell(BspViewCell *viewCell)
{
        mViewCell = viewCell;
}

bool BspLeaf::IsLeaf() const 
{ 
        return true; 
}

void BspLeaf::GenerateViewCell(const BoundedRayContainer &rays, 
                                                           int &sampleContributions,
                                                           int &contributingSamples)
{
        sampleContributions = 0;
        contributingSamples = 0;

        mViewCell = dynamic_cast<BspViewCell *>(ViewCell::Generate());
        mViewCell->mBspLeaves.push_back(this);

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

        // add contributions from samples to the PVS
        for (it = rays.begin(); it != it_end; ++ it)
        {
                int contribution = 0;
                Ray *ray = (*it)->mRay;
                        
                if (!ray->intersections.empty())
                {       
            contribution += mViewCell->GetPvs().AddSample(ray->intersections[0].mObject);
                }

                contribution += mViewCell->GetPvs().AddSample(ray->sourceObject.mObject);

                if (contribution > 0)
                {
                        sampleContributions += contribution;
                        ++ contributingSamples;
                }

                ray->bspLeaves.push_back(this);
        }
}


/****************************************************************/
/*                  class BspTree implementation                */
/****************************************************************/

BspTree::BspTree(ViewCell *viewCell): 
mRootCell(viewCell), mRoot(NULL), mGenerateViewCells(false),
mStorePiercingRays(true)
{
        Randomize(); // initialise random generator for heuristics
}
 
const BspTreeStatistics &BspTree::GetStatistics() const 
{
        return mStat;
}

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

        app << setprecision(4);

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

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

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

        app << "#N_SPLITS ( Number of splits )\n" << splits << "\n";

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

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

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

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

        app << "#N_INPUT_POLYGONS (number of input polygons )\n" <<     polys << endl;

        app << "#N_PVS: " << pvs << endl;

        app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
                 (polys + splits) / (double)polys << endl;
        
        app << "===== END OF BspTree statistics ==========\n";
}

BspTree::~BspTree()
{
        std::stack<BspNode *> tStack;

        tStack.push(mRoot);

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

            tStack.pop();
        
                if (!node->IsLeaf())
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

                        // push the children on the stack (there are always two children)
                        tStack.push(interior->GetBack());
                        tStack.push(interior->GetFront());
                }

                DEL_PTR(node);
        }
}

void BspTree::InsertViewCell(ViewCell *viewCell)
{
        PolygonContainer *polys = new PolygonContainer();

        // extract polygons that guide the split process
        mStat.polys += AddMeshToPolygons(viewCell->GetMesh(), *polys, viewCell);
        mBox.Include(viewCell->GetBox()); // add to BSP aabb

        InsertPolygons(polys);
}

void BspTree::InsertPolygons(PolygonContainer *polys)
{       
        std::stack<BspTraversalData> tStack;
        
        // traverse existing tree or create new tree
    if (!mRoot)
                mRoot = new BspLeaf();

        tStack.push(BspTraversalData(mRoot, polys, 0, mRootCell, new BoundedRayContainer()));

        while (!tStack.empty())
        {
                // filter polygons donw the tree
                BspTraversalData tData = tStack.top();
            tStack.pop();
                        
                if (!tData.mNode->IsLeaf())
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(tData.mNode);

                        //-- filter view cell polygons down the tree until a leaf is reached
                        if (!tData.mPolygons->empty())
                        {
                                PolygonContainer *frontPolys = new PolygonContainer();
                                PolygonContainer *backPolys = new PolygonContainer();
                                PolygonContainer coincident;

                                int splits = 0;
                
                                // split viewcell polygons with respect to split plane
                                splits += interior->SplitPolygons(*tData.mPolygons, 
                                                                                                  *frontPolys, 
                                                                                                  *backPolys,
                                                                                                  coincident,
                                                                                                  sStoreSplitPolys);
                                
                                // extract view cells associated with the split polygons
                                ViewCell *frontViewCell = mRootCell;
                                ViewCell *backViewCell = mRootCell;
        
                                if (!mGenerateViewCells)
                                {
                                        ExtractViewCells(&backViewCell,
                                                                         &frontViewCell,
                                                                         coincident, 
                                                                         interior->mPlane, 
                                                                         frontPolys->empty(), 
                                                                         backPolys->empty());
                                }

                                // don't need coincident polygons anymore
                                interior->ProcessPolygons(&coincident, sStoreSplitPolys);

                                mStat.splits += splits;

                                // push the children on the stack
                                tStack.push(BspTraversalData(interior->GetFront(), 
                                                                                         frontPolys, 
                                                                                         tData.mDepth + 1, 
                                                                                         frontViewCell,
                                                                                         tData.mRays));

                                tStack.push(BspTraversalData(interior->GetBack(), 
                                                                                         backPolys, 
                                                                                         tData.mDepth + 1, 
                                                                                         backViewCell,
                                                                                         new BoundedRayContainer()));
                        }

                        // cleanup
                        DEL_PTR(tData.mPolygons);
                }
                else
                {
                        // reached leaf => subdivide current viewcell
                        BspNode *subRoot = Subdivide(tStack, tData);
                }
        }
}

int BspTree::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())
                {
                        poly->mParent = parent; // set parent intersectable
                        polys.push_back(poly);
                }
                else
                        DEL_PTR(poly);
        }
        return (int)mesh->mFaces.size();
}

int BspTree::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 BspTree::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 BspTree::Construct(const ViewCellContainer &viewCells)
{
        mStat.nodes = 1;
        mBox.Initialize();      // initialise bsp tree bounding box

        // copy view cell meshes into one big polygon soup
        PolygonContainer *polys = new PolygonContainer();
        mStat.polys = AddToPolygonSoup(viewCells, *polys);

        // construct tree from the view cell polygons
        Construct(polys, new BoundedRayContainer());
}


void BspTree::Construct(const ObjectContainer &objects)
{
        mStat.nodes = 1;
        mBox.Initialize();      // initialise bsp tree bounding box
        
        PolygonContainer *polys = new PolygonContainer();
        
        // copy mesh instance polygons into one big polygon soup
        mStat.polys = AddToPolygonSoup(objects, *polys);

        // construct tree from polygon soup
        Construct(polys, new BoundedRayContainer());
}

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

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

        long startTime = GetTime();

        Debug << "**** Extracting polygons from rays ****\n";

        std::map<Face *, Polygon3 *> facePolyMap;

        //-- extract polygons intersected by the rays
        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
        {
                Ray *ray = *rit;
        
                // get ray-face intersection. Store polygon representing the rays together
                // with rays intersecting the face.
                if (!ray->intersections.empty())
                {
                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->intersections[0].mObject);
                        Face *face = obj->GetMesh()->mFaces[ray->intersections[0].mFace];

                        std::map<Face *, Polygon3 *>::iterator it = facePolyMap.find(face);

                        if (it != facePolyMap.end()) 
                        {
                                //store rays if needed for heuristics
                                if (sSplitPlaneStrategy & BLOCKED_RAYS)
                                        (*it).second->mPiercingRays.push_back(ray);
                        } 
                        else
                        {       //store rays if needed for heuristics
                                Polygon3 *poly = new Polygon3(face, obj->GetMesh());
                                poly->mParent = obj;
                                polys->push_back(poly);

                                if (sSplitPlaneStrategy & BLOCKED_RAYS)
                                        poly->mPiercingRays.push_back(ray);

                                facePolyMap[face] = poly;
                        }
                }
        }
        
        facePolyMap.clear();

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

        //-- store rays
        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
        {
                Ray *ray = *rit;
        ray->SetId(-1); // reset id

                float minT, maxT;
                if (BoundRay(*ray, minT, maxT))
                        rays->push_back(new BoundedRay(ray, minT, maxT));
        }

        mStat.polys = (int)polys->size();

        Debug << "**** Finished polygon extraction ****" << endl;
        Debug << (int)polys->size() << " polys extracted from " << (int)rays->size() << " rays" << endl;
        Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;

        Construct(polys, rays);
}

void BspTree::Construct(PolygonContainer *polys, BoundedRayContainer *rays)
{
        std::stack<BspTraversalData> tStack;

        mRoot = new BspLeaf();

        BspTraversalData tData(mRoot, polys, 0, mRootCell, rays);
        tStack.push(tData);

        long startTime = GetTime();
        cout << "**** Contructing bsp tree ****\n";

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

                // subdivide leaf node
                BspNode *subRoot = Subdivide(tStack, tData);
        }

        cout << "**** Finished tree construction ****\n";
        Debug << "BSP tree contruction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
}


BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
{
        //-- terminate traversal  
        if (((int)tData.mPolygons->size() <= sTermMaxPolygons) || 
                ((int)tData.mRays->size() <= sTermMaxRays) ||
                (tData.mDepth >= sTermMaxDepth))
                
        {
                BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
        
                // generate new view cell for each leaf
                if (mGenerateViewCells)
                {
                        int conSamp, sampCon;
                        leaf->GenerateViewCell(*tData.mRays, conSamp, sampCon);
                        
                        mStat.contributingSamples += conSamp;
                        mStat.sampleContributions += sampCon;
                }
                else
                {
                        // add view cell to leaf
                        leaf->SetViewCell(dynamic_cast<BspViewCell *>(tData.mViewCell));
                        leaf->mViewCell->mBspLeaves.push_back(leaf);
                }

                EvaluateLeafStats(tData);
                
                //-- clean up
                
                // remaining polygons are discarded or added to node
                leaf->ProcessPolygons(tData.mPolygons, sStoreSplitPolys);
                DEL_PTR(tData.mPolygons);
                // discard rays
                CLEAR_CONTAINER(*tData.mRays);
                DEL_PTR(tData.mRays);

                return leaf;
        }

        //-- continue subdivision

        PolygonContainer coincident;
        
        PolygonContainer *frontPolys = new PolygonContainer();
        PolygonContainer *backPolys = new PolygonContainer();

        BoundedRayContainer *frontRays = new BoundedRayContainer();
        BoundedRayContainer *backRays = new BoundedRayContainer();
        
        // create new interior node and two leaf nodes
        BspInterior *interior = SubdivideNode(dynamic_cast<BspLeaf *>(tData.mNode), 
                                                                                  *tData.mPolygons,
                                                                                  *frontPolys,
                                                                                  *backPolys,
                                                                                  coincident,
                                                                                  *tData.mRays,
                                                                                  *frontRays,
                                                                                  *backRays);
        ViewCell *frontViewCell = mRootCell;
        ViewCell *backViewCell = mRootCell;

#ifdef _DEBUG   
        if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
        {       for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
                        Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
                Debug << endl;}
#endif

        // extract view cells from coincident polygons according to plane normal
    // only if front or back polygons are empty
        if (!mGenerateViewCells)
        {
                ExtractViewCells(&backViewCell,
                                                 &frontViewCell,
                                                 coincident,
                                                 interior->mPlane,
                                                 backPolys->empty(),
                                                 frontPolys->empty());
        }
        
        // don't need coincident polygons anymore
        interior->ProcessPolygons(&coincident, sStoreSplitPolys);

        // push the children on the stack
        tStack.push(BspTraversalData(interior->GetBack(), 
                                                                 backPolys, 
                                                                 tData.mDepth + 1, 
                                                                 backViewCell, 
                                                                 backRays));

        tStack.push(BspTraversalData(interior->GetFront(), 
                                                                 frontPolys, 
                                                                 tData.mDepth + 1, 
                                                                 frontViewCell, 
                                                                 frontRays));

        // cleanup
        DEL_PTR(tData.mNode);
        DEL_PTR(tData.mPolygons);
        DEL_PTR(tData.mRays);

        return interior;
}

void BspTree::ExtractViewCells(ViewCell **backViewCell, 
                                                           ViewCell **frontViewCell,
                                                           const PolygonContainer &coincident,
                                                           const Plane3 splitPlane, 
                                                           const bool extractBack,
                                                           const bool extractFront) const
{
        bool foundFront = !extractFront, foundBack = !extractBack;

        PolygonContainer::const_iterator it = 
                coincident.begin(), it_end = coincident.end();

        //-- find first view cells in front and back leafs
        for (; !(foundFront && foundBack) && (it != it_end); ++ it)
        {
                if (DotProd((*it)->GetNormal(), splitPlane.mNormal) > 0)
                {
                        *backViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
                        foundBack = true;
                }
                else
                {
                        *frontViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
                        foundFront = true;
                }
        }
        //if (extractFront && foundFront)       Debug << "front view cell: " << *frontViewCell << endl;
        //if (extractBack && foundBack) Debug << "back view cell: " << *backViewCell << endl;
}

BspInterior *BspTree::SubdivideNode(BspLeaf *leaf, 
                                                                        PolygonContainer &polys,
                                                                        PolygonContainer &frontPolys,
                                                                        PolygonContainer &backPolys,
                                                                        PolygonContainer &coincident,
                                                                        BoundedRayContainer &rays,
                                                                        BoundedRayContainer &frontRays,
                                                                        BoundedRayContainer &backRays)
{
        mStat.nodes += 2;
        
        // select subdivision plane
        BspInterior *interior = 
                new BspInterior(SelectPlane(leaf, polys, rays)); 

#ifdef _DEBUG
        Debug << interior << endl;
#endif
        
        // subdivide rays into front and back rays
        SplitRays(interior->mPlane, rays, frontRays, backRays);
        
        // subdivide polygons with plane
        mStat.splits += interior->SplitPolygons(polys, 
                                                    frontPolys, 
                                                                                        backPolys, 
                                                                                        coincident, 
                                                                                        sStoreSplitPolys);

        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));
        
        return interior;
}

void BspTree::SortSplitCandidates(const PolygonContainer &polys, 
                                                                  const int axis, 
                                                                  vector<SortableEntry> &splitCandidates) const
{
  splitCandidates.clear();
  
  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 BspTree::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 + sSplitBorder * (maxBox - minBox);
        float maxBand = minBox + (1.0f - sSplitBorder) * (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;
                        }
                }
        }
  
        float oldCost = (float)polys.size();
        float newCost = sCt_div_ci + minSum / boxArea;
        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 BspTree::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;
                float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
                
                if (r < costRatio)
                {
                        costRatio = r;
                        axis = i;
                        position = p;
                }
        }
        
        if (costRatio >= sMaxCostRatio)
                return false;

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

        return true;
}

Plane3 BspTree::SelectPlane(BspLeaf *leaf, 
                                                        PolygonContainer &polys, 
                                                        const BoundedRayContainer &rays)
{
        if (polys.empty())
        {
                Debug << "Warning: No autopartition polygon candidate available\n";
        
                // return axis aligned split
                AxisAlignedBox3 box;
                box.Initialize();
        
                // create bounding box of region
                Polygon3::IncludeInBox(polys, box);

                const int axis = box.Size().DrivingAxis();
                const Vector3 position = (box.Min()[axis] + box.Max()[axis])*0.5f;

                Vector3 norm(0,0,0); norm[axis] = 1.0f;
                return Plane3(norm, position);
        }
        
        if ((sSplitPlaneStrategy & AXIS_ALIGNED) &&
                ((int)polys.size() > sTermMaxPolysForAxisAligned) &&
                ((int)rays.size() > sTermMaxRaysForAxisAligned) &&
                ((sTermMaxObjectsForAxisAligned < 0) || 
                  (Polygon3::CountIntersectables(polys) > sTermMaxObjectsForAxisAligned)))
        {
                Plane3 plane;
                if (SelectAxisAlignedPlane(plane, polys))
                        return plane;
        }

        // simplest strategy: just take next polygon
        if (sSplitPlaneStrategy & RANDOM_POLYGON)
        {
                return polys[Random((int)polys.size())]->GetSupportingPlane();
        }

        // use heuristics to find appropriate plane
        return SelectPlaneHeuristics(polys, rays, sMaxCandidates);
}

Plane3 BspTree::SelectPlaneHeuristics(PolygonContainer &polys,
                                                                          const BoundedRayContainer &rays,
                                                                          const int maxTests)
{
        float lowestCost = MAX_FLOAT;
        int bestPlaneIdx = 0;
        
        int limit = maxTests > 0 ? Min((int)polys.size(), maxTests) : (int)polys.size();
        
        int candidateIdx = limit;
        
        for (int i = 0; i < limit; ++ i)
        {
                candidateIdx = GetNextCandidateIdx(candidateIdx, polys);

                Plane3 candidatePlane = polys[candidateIdx]->GetSupportingPlane();
                
                // evaluate current candidate
                float candidateCost = SplitPlaneCost(candidatePlane, polys, rays);
                        
                if (candidateCost < lowestCost)
                {
                        bestPlaneIdx = candidateIdx;
                        lowestCost = candidateCost;
                }
        }
        
        //Debug << "Plane lowest cost: " << lowestCost << endl;
        return polys[bestPlaneIdx]->GetSupportingPlane();
}

int BspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
{
        const int candidateIdx = Random(currentIdx --);

        // swap candidates to avoid testing same plane 2 times
        std::swap(polys[currentIdx], polys[candidateIdx]);

        /*Polygon3 *p = polys[candidateIdx];
        polys[candidateIdx] = polys[currentIdx];
        polys[currentIdx] = p;*/
        
        return currentIdx;
        //return Random((int)polys.size());
}

float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
                                                          const PolygonContainer &polys) const
{
        float val = 0;

        float sumBalancedPolys = 0;
        float sumSplits = 0;
        float sumPolyArea = 0;
        float sumBalancedViewCells = 0;
        float sumBlockedRays = 0;
        float totalBlockedRays = 0;
        //float totalArea = 0;
        int totalViewCells = 0;

        // need three unique ids for each type of view cell
        // for balanced view cells criterium
        ViewCell::NewMail();
        const int backId = ViewCell::sMailId;
        ViewCell::NewMail();
        const int frontId = ViewCell::sMailId;
        ViewCell::NewMail();
        const int frontAndBackId = ViewCell::sMailId;

        PolygonContainer::const_iterator it, it_end = polys.end();

        for (it = polys.begin(); it != it_end; ++ it)
        {
                const int classification = (*it)->ClassifyPlane(candidatePlane);

                if (sSplitPlaneStrategy & BALANCED_POLYS)
                        sumBalancedPolys += sBalancedPolysTable[classification];
                
                if (sSplitPlaneStrategy & LEAST_SPLITS)
                        sumSplits += sLeastPolySplitsTable[classification];

                if (sSplitPlaneStrategy & LARGEST_POLY_AREA)
                {
                        if (classification == Polygon3::COINCIDENT)
                                sumPolyArea += (*it)->GetArea();
                        //totalArea += area;
                }
                
                if (sSplitPlaneStrategy & BLOCKED_RAYS)
                {
                        const float blockedRays = (float)(*it)->mPiercingRays.size();
                
                        if (classification == Polygon3::COINCIDENT)
                                sumBlockedRays += blockedRays;
                        
                        totalBlockedRays += blockedRays;
                }

                // assign view cells to back or front according to classificaion
                if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
                {
                        MeshInstance *viewCell = (*it)->mParent;
                
                        // assure that we only count a view cell 
                        // once for the front and once for the back side of the plane
                        if (classification == Polygon3::FRONT_SIDE)
                        {
                                if ((viewCell->mMailbox != frontId) && 
                                        (viewCell->mMailbox != frontAndBackId))
                                {
                                        sumBalancedViewCells += 1.0;

                                        if (viewCell->mMailbox != backId)
                                                viewCell->mMailbox = frontId;
                                        else
                                                viewCell->mMailbox = frontAndBackId;
                                        
                                        ++ totalViewCells;
                                }
                        }
                        else if (classification == Polygon3::BACK_SIDE)
                        {
                                if ((viewCell->mMailbox != backId) &&
                                    (viewCell->mMailbox != frontAndBackId))
                                {
                                        sumBalancedViewCells -= 1.0;

                                        if (viewCell->mMailbox != frontId)
                                                viewCell->mMailbox = backId;
                                        else
                                                viewCell->mMailbox = frontAndBackId; 

                                        ++ totalViewCells;
                                }
                        }
                }
        }

        // all values should be approx. between 0 and 1 so they can be combined
        // and scaled with the factors according to their importance
        if (sSplitPlaneStrategy & BALANCED_POLYS)
                val += sBalancedPolysFactor * fabs(sumBalancedPolys) / (float)polys.size();

        if (sSplitPlaneStrategy & LEAST_SPLITS)    
                val += sLeastSplitsFactor * sumSplits / (float)polys.size();

        if (sSplitPlaneStrategy & LARGEST_POLY_AREA) 
                // HACK: polys.size should be total area so scaling is between 0 and 1
                val += sLargestPolyAreaFactor * (float)polys.size() / sumPolyArea;

        if (sSplitPlaneStrategy & BLOCKED_RAYS)
                if (totalBlockedRays != 0)
                        val += sBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;

        if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
                if (totalViewCells != 0)
                        val += sBalancedViewCellsFactor * fabs(sumBalancedViewCells) / 
                                (float)totalViewCells;
        
        return val;
}

bool BspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
{
        maxT = 1e6;
        minT = 0;
        
        // test with tree bounding box
        if (!mBox.GetMinMaxT(ray, &minT, &maxT))
                return false;

        if (minT < 0) // start ray from origin
                minT = 0;

        // bound ray or line segment
        if ((ray.GetType() == Ray::LOCAL_RAY) && 
            !ray.intersections.empty() && 
                (ray.intersections[0].mT <= maxT))
        {
                maxT = ray.intersections[0].mT;
        }

        return true;
}

float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
                                                          const BoundedRayContainer &rays) const
{
        float val = 0;

        float sumBalancedRays = 0;
        float sumRaySplits = 0;
        float pvsSize = 0;
        float totalSize = 0;

        // need three unique ids for small pvs criterium
        Intersectable::NewMail();
        const int backId = ViewCell::sMailId;
        Intersectable::NewMail();
        const int frontId = ViewCell::sMailId;
        Intersectable::NewMail();
        const int frontAndBackId = ViewCell::sMailId;
        //Debug << "f " << frontId <<  " b " << backId << " fb " << frontAndBackId << endl;

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

        for (rit = rays.begin(); rit != rays.end(); ++ rit)
        {
                Ray *ray = (*rit)->mRay;
                const float minT = (*rit)->mMinT;
                const float maxT = (*rit)->mMaxT;

                const int cf = 
                        ray->ClassifyPlane(candidatePlane, minT, maxT);

                if (sSplitPlaneStrategy & LEAST_RAY_SPLITS)
                {
                        sumBalancedRays += sBalancedRaysTable[cf];
                }
                
                if (sSplitPlaneStrategy & BALANCED_RAYS)
                {
                        sumRaySplits += sLeastRaySplitsTable[cf];
                }

                if (sSplitPlaneStrategy & PVS)
                {
                        vector<Ray::Intersection>::const_iterator it, 
                                it_end = ray->intersections.end();

                        if (!ray->intersections.empty())
                        {

                                // assure that we only count a object 
                                // once for the front and once for the back side of the plane
                                pvsSize += PvsValue(*ray->intersections[0].mObject, 
                                                                        cf, 
                                                                        frontId, 
                                                                        backId, 
                                                                        frontAndBackId);

                                // always add 2 for each object (= maximal large PVS)
                                //if (inc > 0.01)       totalSize += 2.0;
                        }
                        //todo emerging obj
                }
        }

        if (sSplitPlaneStrategy & LEAST_RAY_SPLITS)
                if (!rays.empty())
                        val += sLeastRaySplitsFactor * sumRaySplits / (float)rays.size();

        if (sSplitPlaneStrategy & BALANCED_RAYS)
                if (!rays.empty())
                        val += sBalancedRaysFactor * fabs(sumBalancedRays) / (float)rays.size();

        if (sSplitPlaneStrategy & PVS)
                if (!rays.empty()) // HACK (should be maximal possible pvs)
                        val += sPvsFactor * pvsSize / (float)rays.size(); 

        //Debug << "pvs: " << pvsSize << " val " << val << endl;
        return val;
}


float BspTree::PvsValue(Intersectable &obj, 
                                            const int cf, 
                                            const int frontId, 
                                            const int backId, 
                                            const int frontAndBackId) const
{
        float pvsVal = 0;

        if (cf == Ray::COINCIDENT)
                return pvsVal;

        if (cf == Ray::FRONT)
        {
                if ((obj.mMailbox != frontId) && 
                        (obj.mMailbox != frontAndBackId))
                        pvsVal = 1.0;

                if (obj.mMailbox != backId)
                        obj.mMailbox = frontId;
                else
                        obj.mMailbox = frontAndBackId;                                  
        }
        else if (cf == Ray::BACK)
        {
                if ((obj.mMailbox != backId) &&
                        (obj.mMailbox != frontAndBackId))
                {
                        pvsVal = 1.0;

                        if (obj.mMailbox != frontId)
                                obj.mMailbox = backId;
                        else
                                obj.mMailbox = frontAndBackId; 
                }
        }
        else if ((cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT))
        {
                if ((obj.mMailbox == backId) || (obj.mMailbox == frontId))
                        pvsVal = 1.0;
                else if (obj.mMailbox != frontAndBackId)
                        pvsVal = 2.0;
                        
                obj.mMailbox = frontAndBackId;
        }

        return pvsVal;
}

float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
                                                          const PolygonContainer &polys,                                                          
                                                          const BoundedRayContainer &rays) const
{
        float val = 0;

        if (sSplitPlaneStrategy & VERTICAL_AXIS)
        {
                Vector3 tinyAxis(0,0,0); tinyAxis[mBox.Size().TinyAxis()] = 1.0f;
                // we put a penalty on the dot product between the "tiny" vertical axis
                // and the split plane axis
                val += sVerticalSplitsFactor * 
                           fabs(DotProd(candidatePlane.mNormal, tinyAxis));
        }

        // the following criteria loop over all polygons to find the cost value
        if ((sSplitPlaneStrategy & BALANCED_POLYS)      ||
                (sSplitPlaneStrategy & LEAST_SPLITS)        ||
                (sSplitPlaneStrategy & LARGEST_POLY_AREA)   ||
                (sSplitPlaneStrategy & BALANCED_VIEW_CELLS) ||
                (sSplitPlaneStrategy & BLOCKED_RAYS))
        {
        val += SplitPlaneCost(candidatePlane, polys);
        }

        // the following criteria loop over all rays to find the cost value
        if ((sSplitPlaneStrategy & BALANCED_RAYS)      ||
                (sSplitPlaneStrategy & LEAST_RAY_SPLITS)   ||
                (sSplitPlaneStrategy & PVS))
        {
                val += SplitPlaneCost(candidatePlane, rays);
        }

        // return linear combination of the sums
        return val;
}

void BspTree::ParseEnvironment()
{
        //-- parse bsp cell tree construction method
        char constructionMethodStr[60];
        
        environment->GetStringValue("BspTree.Construction.input", constructionMethodStr);

        sConstructionMethod = FROM_INPUT_VIEW_CELLS;
        
        if (strcmp(constructionMethodStr, "fromViewCells") == 0)
                sConstructionMethod = FROM_INPUT_VIEW_CELLS;
        else if (strcmp(constructionMethodStr, "fromSceneGeometry") == 0)
                sConstructionMethod = FROM_SCENE_GEOMETRY;
        else if (strcmp(constructionMethodStr, "fromRays") == 0)
                sConstructionMethod = FROM_RAYS;
        else 
        {
                cerr << "Wrong construction method " << constructionMethodStr << endl;
                exit(1);
    }

        Debug << "Construction method: " << constructionMethodStr << endl;

        //-- termination criteria for autopartition
        environment->GetIntValue("BspTree.Termination.maxDepth", sTermMaxDepth);
        environment->GetIntValue("BspTree.Termination.maxPolygons", sTermMaxPolygons);
        environment->GetIntValue("BspTree.Termination.maxRays", sTermMaxRays);

        //-- termination criteria for axis aligned split
        environment->GetFloatValue("BspTree.Termination.AxisAligned.ct_div_ci", sCt_div_ci);
        environment->GetFloatValue("BspTree.Termination.AxisAligned.maxCostRatio", sMaxCostRatio);
        environment->GetIntValue("BspTree.Termination.AxisAligned.maxPolys", 
                                                         sTermMaxPolysForAxisAligned);
        environment->GetIntValue("BspTree.Termination.AxisAligned.maxRays", 
                                                         sTermMaxRaysForAxisAligned);
        environment->GetIntValue("BspTree.Termination.AxisAligned.maxObjects", 
                                                         sTermMaxObjectsForAxisAligned);
        //-- partition criteria
        environment->GetIntValue("BspTree.maxCandidates", sMaxCandidates);
        environment->GetIntValue("BspTree.splitPlaneStrategy", sSplitPlaneStrategy);
        environment->GetFloatValue("BspTree.AxisAligned.splitBorder", sSplitBorder);
        
        environment->GetBoolValue("BspTree.storeSplitPolys", sStoreSplitPolys);

        environment->GetFloatValue("BspTree.Construction.sideTolerance", Vector3::sDistTolerance);
        Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;

        environment->GetIntValue("ViewCells.minPvsDif", sMinPvsDif);

    Debug << "BSP max depth: " << sTermMaxDepth << endl;
        Debug << "BSP max polys: " << sTermMaxPolygons << endl;
        Debug << "BSP max candidates: " << sMaxCandidates << endl;
        
        Debug << "Split plane strategy: ";
        if (sSplitPlaneStrategy & RANDOM_POLYGON)
                Debug << "random polygon ";
        if (sSplitPlaneStrategy & AXIS_ALIGNED)
                Debug << "axis aligned ";
        if (sSplitPlaneStrategy & LEAST_SPLITS)     
                Debug << "least splits ";
        if (sSplitPlaneStrategy & BALANCED_POLYS)
                Debug << "balanced polygons ";
        if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
                Debug << "balanced view cells ";
        if (sSplitPlaneStrategy & LARGEST_POLY_AREA)
                Debug << "largest polygon area ";
        if (sSplitPlaneStrategy & VERTICAL_AXIS)
                Debug << "vertical axis ";
        if (sSplitPlaneStrategy & BLOCKED_RAYS)
                Debug << "blocked rays ";
        if (sSplitPlaneStrategy & PVS)
                Debug << "pvs";
        Debug << endl;
}

void BspTree::CollectLeaves(vector<BspLeaf *> &leaves)
{
        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 BspTree::GetBoundingBox() const
{
        return mBox;
}

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

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

        if (data.mDepth >= sTermMaxDepth)
                ++ mStat.maxDepthNodes;
        // store maximal and minimal depth
        if (data.mDepth > mStat.maxDepth)
                mStat.maxDepth = data.mDepth;
        if (data.mDepth < mStat.minDepth)
                mStat.minDepth = data.mDepth;

        // accumulate depth to compute average
        mStat.accumDepth += data.mDepth;

        if (leaf->mViewCell != mRootCell)
        {       
                ++ mStat.viewCells;
                mStat.pvs += leaf->mViewCell->GetPvs().GetSize();
        }
        
#ifdef _DEBUG
        Debug << "BSP stats: "
                  << "Depth: " << data.mDepth << " (max: " << sTermMaxDepth << "), "
                  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << sTermMaxPolygons << "), "
                  << "#rays: " << (int)data.mRays->size() << " (max: " << sTermMaxRays << ")" << endl;
#endif
}

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

        if (!BoundRay(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 = (BspInterior *) node;
                        
                        Plane3 *splitPlane = in->GetPlane();

                        int entSide = splitPlane->Side(entp);
                        int extSide = splitPlane->Side(extp);

                        Vector3 intersection;

                        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 ?
                        {
                                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->mViewCell->Mailed())
                        {
                                ray.bspLeaves.push_back(leaf);
                                leaf->mViewCell->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 BspTree::Export(const string filename)
{
        Exporter *exporter = Exporter::GetExporter(filename);

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

        return false;
}

void BspTree::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)->mViewCell;

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

                        nodeStack.push(interior->mFront);
                        nodeStack.push(interior->mBack);
                }
        }
}

int BspTree::CountViewCellPvs() const
{
        int count = 0;
        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()) 
                {
                        ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->mViewCell;

                        if (!viewCell->Mailed()) 
                        {
                                viewCell->Mail();
                                count += viewCell->GetPvs().GetSize();
                        }
                }
                else 
                {
                        BspInterior *interior = dynamic_cast<BspInterior *>(node);

                        nodeStack.push(interior->mFront);
                        nodeStack.push(interior->mBack);
                }
        }
        return count;
}

bool BspTree::MergeViewCells(BspLeaf *front, BspLeaf *back) const
{
        BspViewCell *viewCell = 
                dynamic_cast<BspViewCell *>(ViewCell::Merge(*front->mViewCell, *back->mViewCell));
        
        if (!viewCell)
                return false;

        // change view cells of all leaves associated with the
        // previous view cells

        BspViewCell *fVc = front->mViewCell;
        BspViewCell *bVc = back->mViewCell;

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

        vector<BspLeaf *>::const_iterator it;
        
        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
        {
                (*it)->SetViewCell(viewCell);
                viewCell->mBspLeaves.push_back(*it);
        }
        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
        {
                (*it)->SetViewCell(viewCell);
                viewCell->mBspLeaves.push_back(*it);
        }
        
        DEL_PTR(fVc);
        DEL_PTR(bVc);

        return true;
}

bool BspTree::ShouldMerge(BspLeaf *front, BspLeaf *back) const
{
        ViewCell *fvc = front->mViewCell;
        ViewCell *bvc = back->mViewCell;

        if ((fvc == mRootCell) || (bvc == mRootCell) || (fvc == bvc))
                return false;

        /*Debug << "merge (" << sMinPvsDif << ")" << endl;
        Debug << "size f: " << fvc->GetPvs().GetSize() << ", "
                  << "size b: " << bvc->GetPvs().GetSize() << endl;
        
        Debug << "diff f: " << fvc->GetPvs().Diff(bvc->GetPvs()) << ", "
                  << "diff b: " << bvc->GetPvs().Diff(fvc->GetPvs()) << endl;*/

    if ((fvc->GetPvs().Diff(bvc->GetPvs()) < sMinPvsDif) &&
            (bvc->GetPvs().Diff(fvc->GetPvs()) < sMinPvsDif))
                return true;

        return false;
}

void BspTree::SetGenerateViewCells(int generateViewCells)
{
        mGenerateViewCells = generateViewCells;
}

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

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

        while (!rays.empty())
        {
                BoundedRay *bRay = rays.back();
                Ray *ray = bRay->mRay;
                rays.pop_back();
                
                const int cf = 
                        ray->ClassifyPlane(plane, bRay->mMinT, bRay->mMaxT);
                
                ray->SetId(cf);

                switch (cf)
                {
                case Ray::COINCIDENT:
                        frontRays.push_back(bRay);
                        break;
                case Ray::BACK:
                        backRays.push_back(bRay);
                        break;
                case Ray::FRONT:
                        frontRays.push_back(bRay);
                        break;
                case Ray::FRONT_BACK:
                        {
                                // find intersection of ray segment with plane
                                const Vector3 extp = ray->Extrap(bRay->mMaxT);
                                const float t = plane.FindT(ray->GetLoc(), extp);
                                
                                const float newT = t * bRay->mMaxT;

                                frontRays.push_back(new BoundedRay(ray, bRay->mMinT, newT));
                                backRays.push_back(new BoundedRay(ray, newT, bRay->mMaxT));

                                DEL_PTR(bRay);

                                ++ splits;
                        }
                        break;
                case Ray::BACK_FRONT:
                        {
                                // find intersection of ray segment with plane
                                const Vector3 extp = ray->Extrap(bRay->mMaxT);
                                const float t = plane.FindT(ray->GetLoc(), extp);
                                const float newT = t * bRay->mMaxT;

                                backRays.push_back(new BoundedRay(ray, bRay->mMinT, newT));
                                frontRays.push_back(new BoundedRay(ray, newT, bRay->mMaxT));
                                
                                DEL_PTR(bRay);

                                ++ splits;
                        }
                        break;
                default:
                        Debug << "Should not come here" << endl;
                        break;
                }
        }

        return splits;
}

void BspTree::ExtractSplitPlanes(BspNode *n, 
                                                                 vector<Plane3 *> &planes, 
                                                                 vector<bool> &sides) 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);

                        planes.push_back(dynamic_cast<BspInterior *>(interior)->GetPlane());
                        sides.push_back(interior->mFront == lastNode);
                }
        }
        while (n);
}

struct Candidate
{
        Polygon3 *mPoly;
        Plane3 mPlane;
        bool mSide;

        Candidate(Polygon3 *poly, const Plane3 &plane, const bool &side):
        mPoly(poly), mPlane(plane), mSide(side)
        {}
};

void BspTree::AddHalfspace(PolygonContainer &cell, 
                                                   vector<Plane3> &planes,
                                                   vector<bool> &sides,
                                                   const Plane3 &halfspace, 
                                                   const bool side) const
{
        Polygon3 *poly = GetBoundingBox().CrossSection(halfspace);

        bool inside = true;

        // polygon is split by all other planes
        for (int i = 0; (i < planes.size()) && inside; ++ i)
        {
                VertexContainer splitPts;
                Polygon3 *frontPoly, *backPoly;

                const int cf = poly->ClassifyPlane(planes[i]);
                        
                // split new polygon with all previous planes
                switch (cf)
                {
                        case Polygon3::SPLIT:
                                frontPoly = new Polygon3();
                                backPoly = new Polygon3();

                                poly->Split(planes[i], *frontPoly, *backPoly, splitPts);
                                DEL_PTR(poly);

                                if(sides[i] == true)
                                {
                                        poly = frontPoly;
                                        DEL_PTR(backPoly);
                                }
                                        else
                                        {
                                                poly = backPoly;
                                                DEL_PTR(frontPoly);
                                        }
                                        inside = true;
                                        break;
                        case Polygon3::BACK_SIDE:
                                if (sides[i]) 
                                        inside = false;
                                break;
                        case Polygon3::FRONT_SIDE:
                                if (!sides[i])
                                        inside = false;
                                break;
                        default:
                                break;
                }
        }

        //-- plane poly splits all other candidates

    // some polygons may fall out => rebuild cell
        vector<Candidate> candidates;

        for (int i = 0; i < (int)cell.size(); ++ i)
        {
                candidates.push_back(Candidate(cell[i], planes[i], sides[i]));
        }

        cell.clear();
        planes.clear();
        sides.clear();

        for (int i = 0; i < (int)candidates.size(); ++ i)
        {
                bool candidateInside = true;

                VertexContainer splitPts;
                Polygon3 *frontPoly, *backPoly;

                const int cf = poly->ClassifyPlane(halfspace);
                        
                // split new polygon with all previous planes
                switch (cf)
                {
                        case Polygon3::SPLIT:
                                frontPoly = new Polygon3();
                                backPoly = new Polygon3();

                                poly->Split(halfspace, *frontPoly, *backPoly, splitPts);
                                DEL_PTR(candidates[i].mPoly);

                                if (sides[i] == true)
                                {
                                        candidates[i].mPoly = frontPoly;
                                        DEL_PTR(backPoly);
                                }
                                else
                                {
                                        candidates[i].mPoly = backPoly;
                                        DEL_PTR(frontPoly);
                                }
                                inside = true;
                                break;
                        case Polygon3::BACK_SIDE:
                                if (candidates[i].mSide) 
                                        candidateInside = false;
                                break;
                        case Polygon3::FRONT_SIDE:
                                if (!candidates[i].mSide)
                                        candidateInside = false;
                                break;
                        default:
                                break;
                }
                
                if (candidateInside)
                {
                        cell.push_back(candidates[i].mPoly);
                        sides.push_back(candidates[i].mSide);
                        planes.push_back(candidates[i].mPlane);
                }
                else
                {
                        DEL_PTR(candidates[i].mPoly);
                }
        }

        //-- finally add polygon
        if (inside)
        {
                cell.push_back(poly);
                planes.push_back(halfspace);
                sides.push_back(side);
        }
        else
        {
                DEL_PTR(poly);
        }
}


void BspTree::ConstructGeometry(BspNode *n, PolygonContainer &cell) const
{
        stack<BspNode *> tStack;

        vector<Plane3 *> planes;
        vector<bool> sides;

        ExtractSplitPlanes(n, planes, sides);

        PolygonContainer candidatePolys;

        // bounded planes are added to the polygons
        for (int i = 0; i < (int)planes.size(); ++ i)
        {
                candidatePolys.push_back(GetBoundingBox().CrossSection(*planes[i]));
        }

        // 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)
        {
                bool inside = true;

                // polygon is split by all other planes
                for (int j = 0; (j < planes.size()) && inside; ++ j)
                {
                        Plane3 *plane = planes[j];

                        if (i == j) // same plane
                                continue;

                        VertexContainer splitPts;
                        Polygon3 *frontPoly, *backPoly;

                        const int cf = candidatePolys[i]->ClassifyPlane(*plane);
                        
                        switch (cf)
                        {
                                case Polygon3::SPLIT:
                                        frontPoly = new Polygon3();
                                        backPoly = new Polygon3();

                                        candidatePolys[i]->Split(*plane, *frontPoly, 
                                                                                         *backPoly, splitPts);
                                        DEL_PTR(candidatePolys[i]);

                                        if(sides[j] == true)
                                        {
                                                candidatePolys[i] = frontPoly;
                                                DEL_PTR(backPoly);
                                        }
                                        else
                                        {
                                                candidatePolys[i] = backPoly;
                                                DEL_PTR(frontPoly);
                                        }
                                        inside = true;
                                        break;

                                case Polygon3::BACK_SIDE:
                                        if (sides[j]) 
                                                inside = false;
                                        break;
                                case Polygon3::FRONT_SIDE:
                                        if (!sides[j])
                                                inside = false;
                                        break;
                                default:
                                        break;
                        }
                }
                
                if (inside)
                        cell.push_back(candidatePolys[i]);
                else
                        DEL_PTR(candidatePolys[i]);
        }
}

void BspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &cell) const
{
        vector<BspLeaf *> leaves = vc->mBspLeaves;

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

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

int BspTree::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 *> planes;
        vector<bool> sides;

        ExtractSplitPlanes(n, planes, sides);

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

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

                                        if ((cf == Polygon3::BACK_SIDE) && sides[i])
                                                isAdjacent = false;
                                        else if ((cf == Polygon3::FRONT_SIDE) && !sides[i])
                                                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->mPlane);

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

BspLeaf *BspTree::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);

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

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

BspLeaf *BspTree::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->mBack);
                        else
                                nodeStack.push(interior->mFront);

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

int BspTree::CountPvs(const BoundedRayContainer &rays) const
{
        int pvsSize = 0;

        BoundedRayContainer::const_iterator rit, rit_end = rays.end();
        vector<BspLeaf *>::const_iterator lit;

        ViewCell::NewMail();

        for (rit = rays.begin(); rit != rit_end; ++ rit)
        {
                Ray *ray = (*rit)->mRay;
                
                for (lit = ray->bspLeaves.begin(); lit != ray->bspLeaves.end(); ++ lit)
                {
                        if ((*lit)->mViewCell->Mailed())
                        {
                                pvsSize += (*lit)->mViewCell->GetPvs().GetSize();
                        }
                }
        }

        return pvsSize;
}