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

#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::sTermMaxRays = -1;

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;

bool BspTree::sStoreSplitPolys = false;

/** Evaluates split plane classification with respect to the plane's 
        contribution for a minimum number splits in the tree.
*/
float BspTree::sLeastSplitsTable[] = {0, 0, 1, 0};
/** Evaluates split plane classification with respect to the plane's 
        contribution for a balanced tree.
*/
float BspTree::sBalancedTreeTable[] = {-1, 1, 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);
}

void BspInterior::SplitRays(RayContainer &rays, 
                                                        RayContainer &frontRays, 
                                                        RayContainer &backRays)
{
        while (!rays.empty())
        {
                Ray *ray = rays.back();
                rays.pop_back();

                //-- ray-plane intersection
                float maxT = 1e6;
                float minT = 0;
        
                // test with tree bounding box
/*              if (!mBox.GetMinMaxT(*ray, &minT, &maxT))
                        continue;
                if (minT < 0) // start ray from origin
                        minT = 0;
*/
                // bound ray
                if ((ray->GetType() == Ray::LOCAL_RAY) && 
                        (!ray->intersections.empty())      &&
                        (ray->intersections[0].mT <= maxT))
                {
                        maxT = ray->intersections[0].mT;
                }
                
                int classification = ray->ClassifyPlane(mPlane, minT, maxT);
                
                switch (classification)
                {
                case Plane3::COINCIDENT:
                        break;
                case Plane3::BACK_SIDE:
                        ray->SetId(BACK_RAY);
                        backRays.push_back(ray);
                        break;
                case Plane3::FRONT_SIDE:
                        ray->SetId(FRONT_RAY);
                        frontRays.push_back(ray);
                        break;
                case Plane3::SPLIT:
                        ray->SetId(SPLIT_RAY);
                        frontRays.push_back(ray);
                        backRays.push_back(ray);
                        //extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &maxt);
                        break;
                default:
                        break;
                }
        }
}

int BspInterior::SplitPolygons(PolygonContainer &polys, 
                                                           PolygonContainer &frontPolys, 
                                                           PolygonContainer &backPolys, 
                                                           PolygonContainer &coincident,
                                                           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
                int classification = poly->ClassifyPlane(mPlane);

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

                switch (classification)
                {
                        case Plane3::COINCIDENT:
                                coincident.push_back(poly);
                                break;                  
                        case Plane3::FRONT_SIDE:        
                                frontPolys.push_back(poly);
                                break;
                        case Plane3::BACK_SIDE:
                                backPolys.push_back(poly);
                                break;
                        case Plane3::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
                                if (poly->mPiercingRays)
                                        InheritRays(*poly, *front_piece, *back_piece);

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

void BspInterior::InheritRays(const Polygon3 &poly, 
                                                      Polygon3 &front_piece, 
                                                      Polygon3 &back_piece)
{
        RayContainer *rays = poly.mPiercingRays;

        RayContainer::const_iterator it, it_end = rays->end();

        for (it = rays->begin(); it != it_end; ++ it)
        {
                switch((*it)->GetId())
                {
                case BACK_RAY:
                        back_piece.GetPiercingRays()->push_back(*it);
                        break;
                case FRONT_RAY:
                        front_piece.GetPiercingRays()->push_back(*it);
                        break;
                }
        }
}

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

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

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

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

ViewCell *BspLeaf::GetViewCell()
{
        return mViewCell;
}

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

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

int BspLeaf::GenerateViewCell(const RayContainer &rays)
{
        int contributingSamples = 0;

        mViewCell = new ViewCell();

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

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

        return contributingSamples;
}


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

BspTree::BspTree(ViewCell *viewCell): 
mRootCell(viewCell), mRoot(NULL), mGenerateViewCells(false)
{
        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 RayContainer()));

        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 RayContainer()));
                        }

                        // 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);
                poly->mParent = parent; // set parent intersectable
                polys.push_back(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 RayContainer());
}


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 RayContainer());
}

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

        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;
                ray->SetId(-1); // reset id
                rays->push_back(ray);

                // 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()) 
                        {
                                (*it).second->AddPiercingRay(ray);
                        } 
                        else
                        {
                                Polygon3 *poly = new Polygon3(face, obj->GetMesh());
                                polys->push_back(poly);
                                poly->AddPiercingRay(ray);

                                facePolyMap[face] = poly;
                        }
                }
        }
        
        Polygon3::IncludeInBox(*polys, mBox);
        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, RayContainer *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)
                        mStat.contributingSamples += leaf->GenerateViewCell(*tData.mRays);
                else
                        // add view cell to leaf
                        leaf->SetViewCell(tData.mViewCell);
                
                EvaluateLeafStats(tData);

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

                return leaf;
        }

        //-- continue subdivision

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

        RayContainer *frontRays = new RayContainer();
        RayContainer *backRays = new RayContainer();
        
        // 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,
                                                                        RayContainer &rays,
                                                                        RayContainer &backRays,
                                                                        RayContainer &frontRays)
{
        mStat.nodes += 2;
        
        // select subdivision plane
        BspInterior *interior = 
                new BspInterior(SelectPlane(leaf, polys, rays)); 

#ifdef _DEBUG
        Debug << interior << endl;
#endif
        
        // partition rays
        interior->SplitRays(rays, frontRays, backRays);
        
        // split polygons with split 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;
}

Plane3 BspTree::SelectPlane(BspLeaf *leaf, 
                                                        PolygonContainer &polys, 
                                                        const RayContainer &rays)
{
        if (polys.size() == 0)
        {
                Debug << "Warning: No split candidate available\n";
                return Plane3();
        }

        if ((sSplitPlaneStrategy & AXIS_ALIGNED) &&
                ((int)polys.size() > sTermMaxPolysForAxisAligned))
        {
                AxisAlignedBox3 box;
                box.Initialize();
        
                // todo: area subdivision
                Polygon3::IncludeInBox(polys, box);
        
                int objectsBack = 0, objectsFront = 0;
                int axis = 0;
                float costRatio = MAX_FLOAT;
                Vector3 position;

                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)
                {
                        Vector3 norm(0,0,0); norm[axis] = 1.0f;
                        return Plane3(norm, position);
                }
        }

        // 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 RayContainer &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)
{
        int candidateIdx = Random(currentIdx --);

        // swap candidates to avoid testing same plane 2 times
        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;

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

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

                if (sSplitPlaneStrategy & BALANCED_POLYS)
                        sumBalancedPolys += sBalancedTreeTable[classification];
                
                if (sSplitPlaneStrategy & LEAST_SPLITS)
                        sumSplits += sLeastSplitsTable[classification];

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

                // assign view cells to back or front according to classificaion
                if (sSplitPlaneStrategy & BALANCED_VIEW_CELLS)
                {
                        MeshInstance *viewCell = (*it)->mParent;
                
                        if (!viewCell->Mailed())
                        {
                                viewCell->Mail();
                                sumBalancedViewCells += sBalancedTreeTable[classification];

                                ++ 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 * sumBlockedRays / totalBlockedRays;

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

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

        float sumBalancedRays = 0;
        float sumRaySplits = 0;

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

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

                float maxT = 1e6;
                float minT = 0;
        
                // test with tree bounding box
                /*if (!mBox.GetMinMaxT(*ray, &minT, &maxT))
                        continue;
                if (minT < 0) // start ray from origin
                        minT = 0;*/

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

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

                if (sSplitPlaneStrategy & LEAST_RAY_SPLITS)
                {
                        sumBalancedRays += sBalancedTreeTable[classification];
                }
                
                if (sSplitPlaneStrategy & BALANCED_RAYS)
                {
                        sumRaySplits += sLeastSplitsTable[classification];
                }
        }

        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();

        return val;
}

float BspTree::SplitPlaneCost(const Plane3 &candidatePlane, 
                                                          const PolygonContainer &polys,                                                          
                                                          const RayContainer &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))
        {
                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.ct_div_ci", sCt_div_ci);
        environment->GetFloatValue("BspTree.Termination.maxCostRatio", sMaxCostRatio);
        environment->GetIntValue("BspTree.Termination.maxPolysForAxisAligned", 
                                                         sTermMaxPolysForAxisAligned);
        
        //-- partition criteria
        environment->GetIntValue("BspTree.maxCandidates", sMaxCandidates);
        environment->GetIntValue("BspTree.splitPlaneStrategy", sSplitPlaneStrategy);
        environment->GetFloatValue("BspTree.splitBorder", sSplitBorder);
        
        environment->GetBoolValue("BspTree.storeSplitPolys", sStoreSplitPolys);

        environment->GetFloatValue("BspTree.Construction.sideTolerance", Polygon3::sSideTolerance);
        Polygon3::sSideToleranceSqrt = Polygon3::sSideTolerance * Polygon3::sSideTolerance;

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

        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 = 1e6;
        float mint = 0;
        
        Intersectable::NewMail();

        // test with tree bounding box
        if (!mBox.GetMinMaxT(ray, &mint, &maxt))
                return 0;

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

        // bound ray
        if ((ray.GetType() == Ray::LOCAL_RAY) && 
                (!ray.intersections.empty())      &&
            (ray.intersections[0].mT <= maxt))
        {
                maxt = ray.intersections[0].mT;
    }
  
        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();

                        float t = 0;
                        bool coplanar = false;
                
                        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
                        extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &maxt);
                        
                } else // reached leaf => intersection with view cell
                {
                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
      
                        if (!leaf->mViewCell->Mailed())
                        {
                                ray.viewCells.push_back(leaf->mViewCell);
                                leaf->mViewCell->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;

                        //Debug << "leaf: new entrypt: " << entp << " new extp: " << extp << endl;

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

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

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