source: GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.cpp @ 2746

#if TOIMPLEMENT

#include <queue>
#include <stack>
#include "Bvh.h"
#include "yare.h"
#include "PerfTimer.h"
#include "Camera.h"
#include "Settings.h"
#include "Context.h"
#include "TriangleBvh.h"
#include "NodeGeometry.h"
#include "Viewer.h"


#include <fstream>
#include <iostream>
#include <iomanip>

#define INVALID_TEST ((unsigned int)-1)
#define MAX_FLOAT 1e20f

#define TYPE_INTERIOR -2
#define TYPE_LEAF -3
#define USE_TIGHTER_BOUNDS_FOR_ALL 0

const static bool useVbos = true;
unsigned int Bvh::sCurrentVboId = -1;

using namespace std;



/*
3 x---------x 2
  |\         \
  | \         \
  |  \         \
  |     4 x---------x 5
  |   |         | 
0 x   |     x 1 |
   \  |         |
    \ |         |
     \|         |
        7 x---------x 6             
*/

static unsigned int sIndices[] =        
{7, // degenerated
 7, 6, 4, 5, 3, 2, 0, 1,
 1, 4, // degenerated
 4, 3, 7, 0, 6, 1, 5, 2,
 2 // degenerated
};


const static int sNumIndicesPerBox = 20;

/*      Order of vertices
        0 = (0, 0, 0)
        1 = (1, 0, 0)
        2 = (1, 1, 0)
        3 = (0, 1, 0)
        4 = (0, 1, 1)
        5 = (1, 1, 1)
        6 = (1, 0, 1)
        7 = (0, 0, 1)
*/


BvhNode::BvhNode(BvhNode *parent): 
mParent(parent),
mAxis(-1), 
mDepth(0), 
mPlaneMask(0),
mPreferredPlane(0),
mVizBox(NULL),
mLastRenderedFrame(-999),
mFirst(-1), 
mLast(-1),
mNumTestNodes(1),
mTestNodesIdx(-1),
mIndicesPtr(-1)
{
}


void BvhNode::ResetVisibility()
{
        mVisibility.Reset();
        mMeasurements.Reset();

        mLastRenderedFrame = -999;
}


BvhNode::~BvhNode() 
{
        if (mVizBox) mVizBox->unref();
}


Box *BvhNode::GetOrCreateVizBox() 
{
        if (!mVizBox) 
        {
                mVizBox = new Box(mBox); 
                mVizBox->ref();
        }

        return mVizBox;
}


void BvhNode::VisibilityInfo::Reset()
{
        mIsVisible = false;
        mIsFullyVisible = false;

        mVisiblePixels = 0;
        mAssumedVisibleFrames = 0;
        mRemainingVisibleFrames = 0;

        mLastVisitedFrame = -1;
        mLastTestedFrame = -1;
        mTurnedVisibleFrame = 0;

        mTimesInvisible = 0;
        mTimesTested = 0;
        mTimesChangedClassification = 0;
        mAvgChangedClassification = 0;
        mIsFrustumCulled = false;

        mLastTestedVisibleFrame = 0;
        mIsNew = true;
}


void BvhNode::SetFullyVisible(const bool visible)
{
        mVisibility.mIsFullyVisible = visible;
}


BvhLeaf::~BvhLeaf()
{
        if (mTriangleBvh) delete mTriangleBvh;
}


/***********************************************************/
/*                class Bvh implementation                 */
/***********************************************************/



Bvh::Bvh(const GeometryVector &geometry, 
                 DistanceSortedRenderAction *const renderer):
mCamera(NULL), 
mFrameId(-1), 
mRoot(NULL), 
mVertices(NULL),
mIndices(NULL),
mUseTighterBoundsOnlyForLeafTests(false),
mRenderer(renderer),
mTestIndices(NULL),
mCurrentIndicesPtr(0),
mCollectTighterBoundsWithMaxLevel(true)
{
        mGeometry = new NodeGeometry*[geometry.size()];
        mGeometrySize = geometry.size();

        BoundingBox sceneBox;

        // compute scene extent
        for (size_t i = 0; i < mGeometrySize; ++ i)
        {
                mGeometry[i] = geometry[i];
                sceneBox.combine(&mGeometry[i]->mBox);

                mBvhStats.mTriangles += mGeometry[i]->mNumTriangles;
        }


        ////////
        //-- create root

        BvhLeaf *leaf = new BvhLeaf(NULL);
        leaf->mDepth = 0;
        leaf->mFirst = 0;
        leaf->mLast = (int)geometry.size() - 1;
        leaf->mBox = sceneBox;

        OUT1("geometry in root: " << leaf->mLast);

        mRoot = leaf;
        mNumNodes = 1;

        // parameters
        mMaxGeometry = Settings::Global()->get_nvocc_bvh_max_objects();
        mMaxTriangles = Settings::Global()->get_nvocc_bvh_max_triangles();
        mMaxDepth = Settings::Global()->get_nvocc_bvh_max_depth();
        mSplitType = Settings::Global()->get_nvocc_bvh_split_type();

        float sceneArea = sceneBox.getSurface();
        float minAreaRatio = Settings::Global()->get_nvocc_bvh_min_area();

        mMinArea = minAreaRatio * sceneArea;

        mMaxDepthForTestingChildren = Settings::Global()->get_nvocc_bvh_max_depth_for_testing_children();
        mUseTighterBoundsOnlyForLeafTests = (Settings::Global()->get_nvocc_use_tighter_bounds() >= 1);
        
        mAreaRatioThreshold = Settings::Global()->get_nvocc_area_ratio_threshold();
        mVolRatioThreshold = Settings::Global()->get_nvocc_vol_ratio_threshold();
}


Bvh::~Bvh() 
{
        if (mVertices) delete []mVertices; 
        if (mIndices) delete [] mIndices;
        if (mTestIndices) delete [] mTestIndices;
        if (mGeometry) delete [] mGeometry;

        if (mRoot) delete mRoot;
}


void Bvh::Construct()
{
        PerfTimer constructTimer;
        constructTimer.Entry();

        OUT1("Info: Constructing BVH");

        mRoot = SubdivideLeaf((BvhLeaf *)mRoot, 0);

        constructTimer.Exit();

        OUT1("Info: Bvh done in " << constructTimer.TotalCount() << " seconds.");
        OUT1("Info: Number Of BVH Nodes = " << mNumNodes);

        // update stats on leaf nodes
        UpdateNumLeaves(mRoot);
        
        mAvgDepth = 1.0f + log((float)GetNumLeaves()) / log(2.0f);

        BvhLeafContainer leaves;
        leaves.reserve(GetNumLeaves());
        CollectLeaves(mRoot, leaves);

        // compute tighter boundaries for the leaves
        PostProcessLeaves(leaves);
        // compute unique ids
        ComputeIds();
        /// pull up some data from the leaves
        UpdateInteriors(mRoot);
        // recompute the boundaries once
        RecomputeBounds();
        // do stats
        ComputeBvhStats();

        PrintBvhStats();
}


float Bvh::EvaluateSahCost(BvhLeaf *leaf, const int axis, const float position)
{
        // count triangles on the left / right
        int left = 0;
        int right = 0;
        BoundingBox leftBox, rightBox;

        const int candidates = std::max(50, (int)(leaf->CountPrimitives() / 20));

        const float finc = leaf->CountPrimitives() / (float)candidates;

        int i = leaf->mFirst;
        float fi = leaf->mFirst + 0.5f;

        BoundingBox box;

        for (; i <= leaf->mLast;) 
        {
                box = mGeometry[i]->GetBoundingBox();

                if (box.getCenter()[axis] < position) 
                {
                        ++ left;
                        // update the box
                        leftBox.combine(&box);
                } 
                else 
                {
                        ++ right;
                        rightBox.combine(&box);
                }
                 
                fi += finc;
                i = (int)fi;
        }

        float bW = 1.0f;
        float leftRatio = left / (float)leaf->CountPrimitives();
        float rightRatio = right / (float)leaf->CountPrimitives();

        float saLeft = 0.0f;
        float saRight = 0.0f;

        // not a valid split
        if (!left || !right)
                return 1e25f;

        saLeft = leftBox.getSurface();
        saRight = rightBox.getSurface();


        return
                saLeft  * ((1.0f - bW) + bW * leftRatio) + 
                saRight * ((1.0f - bW) + bW * rightRatio);
}


float Bvh::SelectPlaneSah(BvhLeaf *leaf, int &axis, float &minCost)
{
        minCost = MAX_FLOAT;
        float bestPos = minCost;
        int bestAxis = 0;

        //  cout<<"Evaluating axis..."<<endl<<flush;

        const int initialPlanes = 3;

        // initiate the costs
        for (axis = 0; axis < 3; ++ axis) 
        {
                const float size = leaf->mBox.getUpper()[axis] - leaf->mBox.getLower()[axis];

                for (int i = 0; i < initialPlanes; ++ i) 
                {
                        const float pos = leaf->mBox.getLower()[axis] + (i + 1) * size / (initialPlanes + 1);
                        const float cost = EvaluateSahCost(leaf, axis, pos);

                        if (cost < minCost) 
                        {
                                minCost = cost;
                                bestPos = pos;
                                bestAxis = axis;
                        }
                }
        }

        axis = bestAxis;

        //  cout<<axis<<endl<<flush;
        const float shrink = 0.5f;

        // now hierarchically refine the initial interval
        float size = shrink * 
                (leaf->mBox.getUpper()[axis] - leaf->mBox.getLower()[axis]) / (initialPlanes + 1);

        const int steps = 4;
        float cost;

        for (int i = 0; i < steps; ++ i) 
        {
                const float left = bestPos - size;
                const float right = bestPos + size;

                cost = EvaluateSahCost(leaf, axis, left);

                if (cost < minCost) 
                {
                        minCost = cost;
                        bestPos = left;
                }

                cost = EvaluateSahCost(leaf, axis, right);

                if (cost < minCost) 
                {
                        minCost = cost;
                        bestPos = right;
                }

                size = shrink * size;
        }

        //OUT1("best axis: " << axis << " " << bestPos << " " << minCost);

        return bestPos;
}


void Bvh::ResizeVisibilityBuffers(const int maxSize)
{
        std::stack<BvhNode *> nodeStack;

        nodeStack.push(mRoot);

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

                static_cast<BvhLeaf *>(node)->mMeasurements.Reset();

                if (!node->IsLeaf()) 
                {
                        BvhInterior *interior = static_cast<BvhInterior *>(node);
                        
                        nodeStack.push(interior->mBack);
                        nodeStack.push(interior->mFront);
                }
        }
}


int Bvh::SortTriangles(BvhLeaf *leaf, 
                                           const int axis, 
                                           const float position)
{
        int i = leaf->mFirst;
        int j = leaf->mLast;

    while (1) 
        {
                while (mGeometry[i]->GetBoundingBox().getCenter()[axis] < position) 
                        ++ i;
        
                while (position < mGeometry[j]->GetBoundingBox().getCenter()[axis])
                        -- j;

                if (i < j) 
                {
                        std::swap(mGeometry[i], mGeometry[j]);
                        
                        ++ i;
                        -- j;
                }
                else
                {
                        break;
                }
        }

        return j;
}


int Bvh::SortTrianglesSpatialMedian(BvhLeaf *leaf, const int axis)
{
        // spatial median
        float x = leaf->mBox.getCenter()[axis];

        return SortTriangles(leaf, axis, x);
}


int Bvh::SortTrianglesObjectMedian(BvhLeaf *leaf, const int axis, float &pos)
{
        // Now distribute the objects into the subnodes
        // Use QuickMedian sort
        int l = leaf->mFirst;
        int r = leaf->mLast;
        int k = (l + r) / 2;

        while (l < r) 
        {
                int i = l;
                int j = r;

                // get some estimation of the pivot
                pos = mGeometry[(l + r) / 2]->GetBoundingBox().getCenter()[axis];

                while (1) 
                {
                        while ((i  <= leaf->mLast) && (mGeometry[i]->GetBoundingBox().getCenter()[axis] < pos))
                                ++ i;

                        while((j >= leaf->mFirst) && (pos < mGeometry[j]->GetBoundingBox().getCenter()[axis]))
                                -- j;

                        if (i <= j) 
                        {
                                std::swap(mGeometry[i], mGeometry[j]);
                                ++ i;
                                -- j;
                        }
                        else
                        {
                                break;
                        }
                }

                // now check the extents
                if (i >= k)
                        r = j;
                else
                        l = i;
        }

        return k;
}


int Bvh::SortTrianglesSurfaceArea(BvhLeaf *leaf, const float sa)
{
        int i = leaf->mFirst;
        int j = leaf->mLast;

        while(1) 
        {
                while ((i <= j) && (mGeometry[i]->GetBoundingBox().getSurface() < sa)) 
                        ++ i;

                while ((i <= j) && (sa < mGeometry[j]->GetBoundingBox().getSurface()))
                        -- j;

                if (i < j) 
                {
                        swap(mGeometry[i], mGeometry[j]);
                        ++ i;
                        -- j;
                }
                else
                        break;
        }

        return j;
}


bool Bvh::TerminationCriteriaMet(BvhLeaf *leaf) const
{
        const bool terminationCriteriaMet = 
                (leaf->mBox.getSurface() < mMinArea) ||
                (leaf->CountPrimitives() <= mMaxGeometry) ||
                (CountTriangles(leaf) <= mMaxTriangles) ||
                (leaf->mDepth > mMaxDepth);

        return terminationCriteriaMet;
}


BvhNode *Bvh::SubdivideLeaf(BvhLeaf *leaf, const int parentAxis)
{
        if (TerminationCriteriaMet(leaf))
                return leaf;
        
        //int axis = leaf->mBox.MajorAxis();
        int axis = (parentAxis + 1) % 3;

        // position of the split in the partailly sorted array of triangles
        // corresponding to this leaf
        int split = -1;
        const int scale = 20;

        float pos = -1.0f;

        // Spatial median subdivision
        switch (mSplitType) 
        {
        case SPATIAL_MEDIAN:
                split = SortTrianglesSpatialMedian(leaf, axis);

                if ( 
                        ((split - leaf->mFirst) < leaf->CountPrimitives() / scale) ||
                        ((leaf->mLast - split) < leaf->CountPrimitives() / scale) ) 
                {
                        split = SortTrianglesObjectMedian(leaf, axis, pos);
                }

                pos = leaf->mBox.getCenter()[axis];
                break;
        case OBJECT_MEDIAN:
                // Object median subdivision: approximately the same number
                // of objects on the left of the the splitting point.
                split = SortTrianglesObjectMedian(leaf, axis, pos);
                break;
        case SAH:
                {
                        float cost;
                        pos = SelectPlaneSah(leaf, axis, cost);

                        if (pos != MAX_FLOAT)
                        {
                                split = SortTriangles(leaf, axis, pos);
                        }

                        if ((pos == MAX_FLOAT) || (split == leaf->mLast))
                        {
                                split = -1;
                                split = SortTrianglesObjectMedian(leaf, axis, pos);
                        }
                }
                break;
        case SAH_OR_SIZE: 
                {
                        // split by size instead
                        const float saThreshold = 0.2f * leaf->GetBox().getSurface();
                        split = SortTrianglesSurfaceArea(leaf, saThreshold);

                        if ((split == leaf->mLast) || (split == leaf->mFirst - 1)) 
                        {
                                // use SAH
                                float cost;
                                pos = SelectPlaneSah(leaf, axis, cost);

                                if (pos != MAX_FLOAT)
                                        split = SortTriangles(leaf, axis, pos);
                                else
                                        split = SortTrianglesObjectMedian(leaf, axis, pos);
                        }
                        else
                        {
                                // note: no position is computed!!
                                //OUT1("sorted by size");
                        }
                }
                break;
        default:
                OUT1("should not come here");
                break;
        }

        if (1 && ((split == leaf->mLast)))
        {
                // no reduction: we should never come here
                OUT1("error: no reduction " << leaf->CountPrimitives() << " " << leaf->mFirst << " " << split << " " << leaf->mLast);
                return leaf;
        }

        // create two more nodes
        mNumNodes += 2;
        BvhInterior *parent = new BvhInterior(leaf->GetParent());
        BvhLeaf *front = new BvhLeaf(parent);
        
        parent->mAxis = axis;
        parent->mBox = leaf->mBox;
        parent->mDepth = leaf->mDepth;

        parent->mBack = leaf;
        parent->mFront = front;
        parent->mPosition = pos;
        
        // now assign the triangles to the subnodes
        front->mFirst = split + 1;
        front->mLast = leaf->mLast;
        front->mDepth = leaf->mDepth + 1;
        leaf->mLast = split;
        leaf->mDepth = front->mDepth;
        leaf->mParent = parent;
        
        // reset box
        leaf->mBox = BoundingBox();

        UpdateLeafBox(static_cast<BvhLeaf *>(parent->mBack));
        UpdateLeafBox(static_cast<BvhLeaf *>(parent->mFront));

        // some output
        const int n = 500;
        if ((GetNumLeaves() % n) == n - 1)
                OUT1("created " << GetNumLeaves() << " leaves");
        
#if VERBOSE_OUTPUT
        OUT1("box: " << parent->mBox.getUpper() << " " << parent->mBox.getLower());
        OUT1("depth: " << (int)parent->mDepth);
        OUT1("axis: " << axis);
        OUT1("bc="<<((BvhLeaf *)parent->mBack)->Count());
        OUT1("fc="<<((BvhLeaf *)parent->mFront)->Count());
        OUT1("back: " << parent->mBack->mBox.getSurface());
        OUT1("front: " << parent->mFront->mBox.getSurface());
#endif

        // recursively continue subdivision
        parent->mBack = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mBack), axis);
        parent->mFront = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mFront), axis);

        return parent;
}


void Bvh::UpdateLeafBox(BvhLeaf *leaf)
{
        for (int i = leaf->mFirst; i <= leaf->mLast; ++ i) 
        {
                leaf->mBox.combine(&mGeometry[i]->mBox);
        } // for
}


void Bvh::UpdateNumLeaves(BvhNode *node) const
{
        if (node->IsLeaf())
        {
                node->mNumLeaves = 1;
        }
        else
        {
                BvhNode *f = static_cast<BvhInterior *>(node)->GetFront();
                BvhNode *b = static_cast<BvhInterior *>(node)->GetBack();

                UpdateNumLeaves(f);
                UpdateNumLeaves(b);

                node->mNumLeaves = f->mNumLeaves + b->mNumLeaves;
        }
}


void Bvh::UpdateBoxes(BvhNode *node)
{
        if (node->IsLeaf()) 
        {
                UpdateLeafBox(static_cast<BvhLeaf *>(node));
        }
        else 
        {
                BvhInterior *interior = static_cast<BvhInterior *>(node);
        
                UpdateBoxes(interior->mBack);
                UpdateBoxes(interior->mFront);

                node->mBox = interior->mBack->mBox;
                node->mBox.combine(&interior->mFront->mBox);
        }
}


void Bvh::UpdateFullVisibility(BvhNode *node) const
{               
        // node not visited in this frame => no change
        if (node->GetLastVisitedFrame() != mFrameId) 
                return;

        // leaf node: terminate recursion
        if (node->IsLeaf())
        {
                node->SetFullyVisible(node->IsVisible());       
                return;
        }

        BvhInterior *interior = static_cast<BvhInterior *>(node);
        
        // recursive traversal
        UpdateFullVisibility(interior->mBack);
        UpdateFullVisibility(interior->mFront);
        
        interior->SetFullyVisible(interior->mBack->IsFullyVisible() && 
                                                          interior->mFront->IsFullyVisible());
}


void Bvh::PullUpLastVisited(BvhNode *node, const int frameId) const
{               
        BvhNode *parent = node->GetParent();

        while (parent && (parent->GetLastVisitedFrame() != frameId))
        {
                parent->SetLastVisitedFrame(frameId);
                parent = parent->GetParent();
        }
}


void Bvh::MakeParentsVisible(BvhNode *node)
{
        BvhNode *parent = node->GetParent();

        while (parent && (!parent->IsVisible()))
        {
                parent->SetVisible(true);
                parent = parent->GetParent();
        }
}


void Bvh::CollectLeaves(BvhNode *node, BvhLeafContainer &leaves)
{
        stack<BvhNode *> tStack;
        tStack.push(node);

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

                tStack.pop();

                if (!node->IsLeaf())
                {
                        BvhInterior *interior = static_cast<BvhInterior *>(node);

                        tStack.push(interior->mFront);
                        tStack.push(interior->mBack);
                }
                else
                {
                        leaves.push_back(static_cast<BvhLeaf *>(node));
                }
        }
}


void Bvh::CollectNodes(BvhNode *node, BvhNodeContainer &nodes)
{
        stack<BvhNode *> tStack;

        tStack.push(node);

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

                nodes.push_back(node);

                if (!node->IsLeaf())
                {
                        BvhInterior *interior = static_cast<BvhInterior *>(node);

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


typedef pair<BvhNode *, int> tPair;

void Bvh::CollectNodes(BvhNode *root, 
                                           const int depth, 
                                           HierarchyNodeContainer &nodes)
{
        stack<tPair> tStack;
        tStack.push(tPair(root, 0));

        while (!tStack.empty())
        {
                BvhNode *node = tStack.top().first;
                const int d = tStack.top().second;

                tStack.pop();

                // found depth => take this node
                if (d == depth)
                {
                        nodes.push_back(node);
                }
                else if (node->IsLeaf())
                {
                        // found leaf
                        BvhLeaf *leaf = static_cast<BvhLeaf *>(node);

                        // there is no further subdivision on triangle level
                        if (leaf->mTriangleBvh->GetNumNodes() == 1)
                        {
                                nodes.push_back(leaf);
                        }
                        else // more than a root => search in local bvh
                        {
                                leaf->mTriangleBvh->
                                        CollectNodes(leaf->mTriangleBvh->GetRoot(), depth - d, nodes);
                        }
                }
                else // interior
                {
                        BvhInterior *interior = static_cast<BvhInterior *>(node);

                        tStack.push(tPair(interior->mFront, d + 1));
                        tStack.push(tPair(interior->mBack, d + 1));
                }
        }
}


void Bvh::CollectAllNodes(BvhNode *node, HierarchyNodeContainer &nodes)
{
        BvhNodeContainer bvhNodes;

        // collect nodes of geometry bvh
        CollectNodes(mRoot, bvhNodes);

        
        // first collect nodes of object bvh
        BvhNodeContainer::const_iterator lit, lit_end = bvhNodes.end();

        for (lit = bvhNodes.begin(); lit != lit_end; ++ lit)
        {
                nodes.push_back(*lit);
        }

        // collect nodes of local bvh: all except roots, 
        // because they are equivalent to geometry bvh leaves
        for (lit = bvhNodes.begin(); lit != lit_end; ++ lit)
        {
                if ((*lit)->IsLeaf())
                {
                        BvhLeaf *leaf = static_cast<BvhLeaf *>(*lit);

                        TriangleBvhNodeContainer hnodes;
                        leaf->mTriangleBvh->CollectNodes(leaf->mTriangleBvh->GetRoot(), hnodes);

                        TriangleBvhNodeContainer::const_iterator hit, hit_end = hnodes.end();

                        for (hit = hnodes.begin(); hit != hit_end; ++ hit)
                        {
                                TriangleBvhNode *n = *hit;

                                // don't include root of local bvh - it has the same bounds as the leaves
                                if (n != leaf->mTriangleBvh->GetRoot())
                                        nodes.push_back(n);
                        }
                }
        }
}


void Bvh::CollectVisibleLeaves(BvhNode *node, BvhLeafContainer &leaves)
{
        stack<BvhNode *> tStack;
        tStack.push(node);

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

                if (IsWithinViewFrustum(node))
                {
                        if (!node->IsLeaf())
                        {
                                BvhInterior *interior = static_cast<BvhInterior *>(node);

                                tStack.push(interior->mFront);
                                tStack.push(interior->mBack);
                        }
                        else
                        {
                                leaves.push_back(static_cast<BvhLeaf *>(node));
                        }
                }
        }
}


BvhLeaf *Bvh::GetRandomLeaf(BvhNode *node)
{
        stack<BvhNode *> nodeStack;
  
        nodeStack.push(node);

        int mask = rand();

    while (!nodeStack.empty()) 
        {
                BvhNode *node = nodeStack.top();
                nodeStack.pop();
                
                if (node->IsLeaf()) 
                {
                        // random leaf
                        return static_cast<BvhLeaf *>(node);
                } 
                else 
                {
                        BvhInterior *interior = static_cast<BvhInterior *>(node);
                        
                        BvhNode *next;
                        
                        // random decision
                        if (mask & 1)
                                next = interior->mBack;
                        else
                                next = interior->mFront;

                        mask = mask >> 1;

                        nodeStack.push(next);
                }
        }
  
        // should never come here
        return NULL;
}


BvhLeaf *Bvh::GetRandomVisibleLeaf(BvhNode *node)
{
        BvhLeafContainer leaves;
        leaves.reserve(node->GetNumLeaves());

        CollectVisibleLeaves(node, leaves);

        if (leaves.empty())
                return NULL;

        const int r = (int)(((float)rand() / RAND_MAX) * ((float)leaves.size() - 0.5f));

        return leaves[r];
}


void Bvh::CollectGeometry(BvhNode *node, GeometryVector &geometry)
{
        //geometry.reserve(node->CountPrimitives());

        for (int i = node->mFirst; i <= node->mLast; ++ i)
        {
                geometry.push_back(mGeometry[i]);
        }
}


NodeGeometry **Bvh::GetGeometry(BvhNode *node, int &geometrySize)
{
        geometrySize = node->CountPrimitives();
        return mGeometry + node->mFirst;
}


void Bvh::UpdateInteriors(BvhNode *node)
{
        if (!node->IsLeaf())
        {
                BvhInterior *interior = static_cast<BvhInterior *> (node);
                
                // update the indices of the geometry so we can render interiors as well
                UpdateInteriors(interior->GetBack());
                UpdateInteriors(interior->GetFront());

                // update area
                interior->mFirst = min(interior->GetBack()->mFirst, interior->GetFront()->mFirst);
                interior->mLast = max(interior->GetBack()->mLast, interior->GetFront()->mLast);
        
                interior->mArea = interior->GetBack()->mArea + interior->GetFront()->mArea;
        }
}


int     Bvh::IsWithinViewFrustumLocal(BvhNode *node)
{
        bool bIntersect = false;

        if (node->GetParent())
                node->mPlaneMask = node->GetParent()->mPlaneMask;


        ////////
        //-- do the view frustum culling for the planes [mPreferredPlane - 5]

        for (int i = node->mPreferredPlane; i < 6; ++ i)
        {
                //-- do the test only if necessary
                if (node->mPlaneMask & (1 << i))
                {
                        //-- test the n-vertex
                        if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
                        {
                                //-- outside
                                node->mPreferredPlane = i;
                                return 0;
                        }

                        //-- test the p-vertex
                        if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f)
                        {
                                //-- completely inside: children don't need to check against this plane no more
                                node->mPlaneMask^=      1 << i;
                        }
                        else 
                        {
                                bIntersect = true;
                        }
                }
        }

        //////////
        //-- do the view frustum culling for the planes [0 - m_iPreferredPlane)

        for (int i = 0; i < node->mPreferredPlane; ++ i)
        {
                // do the test only if necessary
                if (node->mPlaneMask & (1 << i))
                {
                        //-- test the n-vertex
                        if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
                        {
                                // outside
                                node->mPreferredPlane = i;
                                return 0;
                        }

                        //-- test the p-vertex
                        if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f)
                        {
                                // completely inside: children don't need to check against this plane no more
                                node->mPlaneMask^= 1 << i;
                        }
                        else 
                        {
                                bIntersect = true;
                        }
                }
        }

        return bIntersect ? -1 : 1;
}


int     Bvh::IsWithinViewFrustum(const BoundingBox &box, const int planeMask, const int preferredPlane)
{
        bool bIntersect = false;

        ////////
        //-- do the view frustum culling for the planes [mPreferredPlane - 5]

        for (int i = preferredPlane; i < 6; ++ i)
        {
                //-- do the test only if necessary
                if (planeMask & (1 << i))
                {
                        //-- test the n-vertex
                        if (box.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
                                return 0;

                        //-- test the p-vertex
                        if (!(box.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f))
                                bIntersect = true;
                }
        }

        //////////
        //-- do the view frustum culling for the planes [0 - m_iPreferredPlane)

        for (int i = 0; i < preferredPlane; ++ i)
        {
                // do the test only if necessary
                if (planeMask & (1 << i))
                {
                        //-- test the n-vertex
                        if (box.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
                                // outside
                                return 0;

                        //-- test the p-vertex
                        if (!(box.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f))
                                bIntersect = true;
                }
        }

        return bIntersect ? -1 : 1;
}


int     Bvh::IsWithinViewFrustum(BvhNode *node)
{
        if (node->mCache.mLastFrustumTestedFrameId == mFrameId)
        {
                return node->mCache.mFrustumCulled;
        }

        node->mCache.mLastFrustumTestedFrameId = mFrameId;

        timeViewFrustumCulling.Entry();
        
        if (Settings::Global()->get_nvocc_use_tighter_bounds_for_frustum_culling())
        {
                const int intersect = IsWithinViewFrustumLocal(node);

                if ((node->mNumTestNodes == 1) || (intersect != -1))
                {
                        //OUT1("x");
                        node->mCache.mFrustumCulled = intersect;
                }
                else 
                {
                        // maybe intersecting one of the tighter boxes
                        node->mCache.mFrustumCulled = 0;

                        for (int i = 0; i < node->mNumTestNodes; ++ i)
                        {
                                RenderableHierarchyNode *n = mTestNodes[node->mTestNodesIdx + i];

                                if (IsWithinViewFrustum(n->GetBox(), 
                                                                                node->mPlaneMask, 
                                                                                node->mPreferredPlane))
                                {
                                        node->mCache.mFrustumCulled = -1;
                                        break;
                                }
                        }
                }
        }
        else
        {
                //OUT1("y");
                node->mCache.mFrustumCulled = IsWithinViewFrustumLocal(node);
        }

        timeViewFrustumCulling.Exit();

        return node->mCache.mFrustumCulled;
}


void Bvh::InitFrame(Camera *camera, const int currentFrameId, Viewer *viewer)
{
        // = 0011 1111 which means that at the beginning, all six planes have to frustum culled
        mRoot->mPlaneMask = 0x3f;
        mCamera = camera;

        mFrameId = currentFrameId;

        // to begin with, we must grab the plane equations of the six clipplanes of the viewfrustum
        Matrix4f matViewing, matProjectionView;
        Transform3D     transform;

        camera->GetTransform(transform);
        transform.get(matViewing);

        camera->GetProjectionMatrix(matProjectionView);
        matProjectionView *= matViewing;

        Vec3f vec;
        float fInvLength;


        //////////
        //-- extract the plane equations

        for (int i = 0; i < 4; ++ i)
        {
                mFrustum.mClipPlane[0][i] = matProjectionView[i][3] - matProjectionView[i][0];  // right plane
                mFrustum.mClipPlane[1][i] = matProjectionView[i][3] + matProjectionView[i][0];  // left plane
                mFrustum.mClipPlane[2][i] = matProjectionView[i][3] + matProjectionView[i][1];  // bottom plane
                mFrustum.mClipPlane[3][i] = matProjectionView[i][3] - matProjectionView[i][1];  // top plane
                mFrustum.mClipPlane[4][i] = matProjectionView[i][3] - matProjectionView[i][2];  // far plane
                mFrustum.mClipPlane[5][i] = matProjectionView[i][3] + matProjectionView[i][2];  // near plane
        }

        ////////////
        //-- normalize the coefficients and find the indices of the n- and p-vertices

        for (int i = 0; i < 6; ++ i)
        {
                // the clipping planes look outward the frustum,
                // so distances > 0 mean that a point is outside
                fInvLength      = -1.0f / sqrt( mFrustum.mClipPlane[i][0] * mFrustum.mClipPlane[i][0] + 
                                                                        mFrustum.mClipPlane[i][1] * mFrustum.mClipPlane[i][1] + 
                                                                        mFrustum.mClipPlane[i][2] * mFrustum.mClipPlane[i][2]);

                mFrustum.mClipPlane[i][0] *= fInvLength;
                mFrustum.mClipPlane[i][1] *= fInvLength;
                mFrustum.mClipPlane[i][2] *= fInvLength;
                mFrustum.mClipPlane[i][3] *= fInvLength;

                vec.x = mFrustum.mClipPlane[i][0];
                vec.y = mFrustum.mClipPlane[i][1];
                vec.z = mFrustum.mClipPlane[i][2];

                // n-vertex
                mClipPlaneAABBVertexIndices[i][0] = BoundingBox::getIndexNearestVertex(vec);
                // p-vertex
                mClipPlaneAABBVertexIndices[i][1] = BoundingBox::getIndexFurthestVertex(vec);   
        }


        const float aspect = mCamera->GetAspect();
        const float fov = mCamera->GetFov();
        
        mScale = sqrt(aspect) * 2.0f * tan(fov * mypi / 180.0f);

        mNearPlane = mCamera->GetRealViewDirection();
        mNearPlaneD = - mNearPlane * (mCamera->GetPosition() - Origin3f);
}


float Bvh::CalcDistance(BvhNode *node) const
{
        timeDistance.Entry();

        float distance;

#if USE_TIGHTER_BOUNDS_FOR_ALL
        float minDist = 1e25f;

        //OUT1("testing " << mNumTestNodes << " nodes");
        for (int i = 0; i < node->mNumTestNodes; ++ i)
        {
                RenderableHierarchyNode *hnode = mTestNodes[node->mTestNodesIdx + i];

                if(hnode->mCache.mLastDistanceTestedFrameId < mFrameId)
                {
                        hnode->mCache.mDistance = 
                                hnode->GetBox().getMinVisibleDistance(mNearPlane, mNearPlaneD);
                        hnode->mCache.mLastDistanceTestedFrameId = mFrameId;
                }

                if (hnode->mCache.mDistance < minDist)
                        minDist = hnode->mCache.mDistance;
        }

        distance = minDist;

#else

        if(node->mCache.mLastDistanceTestedFrameId != mFrameId)
        {
                //OUT1("z " << node->mLastDistanceTestedFrameId  << " " << mFrameId);
                node->mCache.mDistance = node->GetBox().getMinVisibleDistance(mNearPlane, mNearPlaneD);
                node->mCache.mLastDistanceTestedFrameId = mFrameId;
        }

        distance = node->mCache.mDistance;
#endif

        timeDistance.Exit();

        return distance;
}


float Bvh::GetSquareDistance(BvhNode *node) const
{
        timeDistance.Entry();

        float distance;

#if USE_TIGHTER_BOUNDS_FOR_ALL
        const Vector3f nearplane = mCamera->GetRealViewDirection();
        const float nearplaneD = - nearplane * (mCamera->GetPosition() - Origin3f);

        float minDist = 1e25f;
        
        //OUT1("testing " << mNumTestNodes << " nodes");
        for (int i = 0; i < node->mNumTestNodes; ++ i)
        {
                RenderableHierarchyNode *hnode = mTestNodes[node->mTestNodesIdx + i];
                const float dist = GetMinSquareDistance(hnode->GetBox());

                if (dist < minDist)
                        minDist = dist;
        }

        distance = minDist;
#else
        distance = GetMinSquareDistance(node->GetBox());
#endif

        timeDistance.Exit();

        return distance;
}


float Bvh::GetMinSquareDistance(const BoundingBox &box) const
{
        float minDist = 1e25f;

        const Point3f *pts = (const Point3f *)box.getVertexData();
        const Point3f pos = mCamera->GetPosition();

        for (int i = 0; i < 8; ++ i)
        {
                const float newDist = DistanceSquared(pts[i], pos);

                if (newDist < minDist)
                        minDist = newDist;
        }

        return minDist;
}


bool Bvh::ExportBinTree(const string &filename)
{
        ofstream stream(filename.c_str(), ifstream::binary);
        
        if (!stream.is_open()) return false;

        OUT1("exporting bvh");

        std::queue<BvhNode *> tStack;
        tStack.push(mRoot);

        while(!tStack.empty())
        {
                BvhNode *node = tStack.front();
                tStack.pop();
                
                if (node->IsLeaf())
                {
                        //OUT1("l");
                        ExportBinLeaf(stream, static_cast<BvhLeaf *>(node));
                }
                else
                {
                        //OUT1("i");
                        BvhInterior *interior = static_cast<BvhInterior *>(node);

                        ExportBinInterior(stream, interior);
                        
                        tStack.push(interior->mFront);
                        tStack.push(interior->mBack);
                }
        }

        OUT1("... finished");

        return true;
}


void Bvh::ExportBinLeaf(ofstream &stream, BvhLeaf *leaf)
{
        int type = TYPE_LEAF;
        int first = leaf->mFirst;
        int last = leaf->mLast;

        stream.write(reinterpret_cast<char *>(&type), sizeof(int));
        stream.write(reinterpret_cast<char *>(&first), sizeof(int));
        stream.write(reinterpret_cast<char *>(&last), sizeof(int));
}


void Bvh::ExportBinInterior(ofstream &stream, BvhInterior *interior)
{
        int type = TYPE_INTERIOR;
        stream.write(reinterpret_cast<char *>(&type), sizeof(int));
        stream.write(reinterpret_cast<char *>(&interior->mAxis), sizeof(char));
        stream.write(reinterpret_cast<char *>(&interior->mPosition), sizeof(float));
}


BvhLeaf *Bvh::ImportBinLeaf(ifstream &stream, BvhInterior *parent)
{
        BvhLeaf *leaf = new BvhLeaf(parent);

        int first, last;

        stream.read(reinterpret_cast<char *>(&first), sizeof(int));
        stream.read(reinterpret_cast<char *>(&last), sizeof(int));

        leaf->mFirst = first;
        leaf->mLast = last;

        return leaf;
}


BvhInterior *Bvh::ImportBinInterior(ifstream &stream, BvhInterior *parent)
{
        BvhInterior *interior = new BvhInterior(parent);

        stream.read(reinterpret_cast<char *>(&interior->mAxis), sizeof(char));
        stream.read(reinterpret_cast<char *>(&interior->mPosition), sizeof(float));

        return interior;
}


BvhNode *Bvh::LoadNextNode(ifstream &stream, BvhInterior *parent)
{
        int nodeType;
        stream.read(reinterpret_cast<char *>(&nodeType), sizeof(int));

        if (nodeType == TYPE_LEAF)
        {
                //OUT1("l");
                return ImportBinLeaf(stream, static_cast<BvhInterior *>(parent));
        }

        if (nodeType == TYPE_INTERIOR)
        {
                //OUT1("i");
                return ImportBinInterior(stream, static_cast<BvhInterior *>(parent));
        }

        return NULL;
}


Bvh *Bvh::LoadFromFile(const string &filename, const GeometryVector &geom)
{
        // export binary version of mesh
        queue<BvhNode *> tStack;
        ifstream stream(filename.c_str(), ifstream::binary);

        if (!stream.is_open()) return NULL;

        OUT1("loading bvh");
        Bvh *bvh = new Bvh();

        bvh->mGeometrySize = geom.size();
        bvh->mGeometry = new NodeGeometry*[bvh->mGeometrySize];

        for (size_t i = 0; i < bvh->mGeometrySize; ++ i)
                bvh->mGeometry[i] = geom[i];

        bvh->mRoot = bvh->LoadNextNode(stream, NULL);

        tStack.push(bvh->mRoot);
        bvh->mNumNodes = 1;

        while(!tStack.empty())
        {
                BvhNode *node = tStack.front();
                tStack.pop();

                if (!node->IsLeaf())
                {
                        bvh->mNumNodes += 2;

                        BvhInterior *interior = static_cast<BvhInterior *>(node);

            BvhNode *front = bvh->LoadNextNode(stream, interior);
                        BvhNode *back = bvh->LoadNextNode(stream, interior);
        
                        interior->mFront = front;
                        interior->mBack = back;

                        front->mDepth = interior->mDepth + 1;
                        back->mDepth = interior->mDepth + 1;

                        tStack.push(front);                     
                        tStack.push(back);
                }
        }

        OUT1("... finished loading " << bvh->mNumNodes << " nodes, updating boxes");

        //adjust bounding boxes
        bvh->UpdateBoxes(bvh->mRoot);
        bvh->UpdateNumLeaves(bvh->mRoot);

        // compute unique ids
        bvh->ComputeIds();
        // update the indices of the geometry so we can render interiors as well
        bvh->UpdateInteriors(bvh->mRoot);
        // do this once so at least the current node bounding box is tested
        bvh->RecomputeBounds();

        return bvh;
}


float Bvh::ComputeScreenSpaceProjection(BvhNode *node) const
{
#if 0
        int projection = 0;

        for (int i = 0; i < node->mNumTestNodes; ++ i)
        {
                RenderableHierarchyNode *n = mTestNodes[node->mTestNodesIdx + i];

                if(hnode->.mCache.mLastProjectionFrameId < mFrameId)
                {
                        hnode->mCache.mProjection = 
                                ComputeScreenSpaceProjection(hnode->GetBox());
                        hnode->mCache.mLastProjectionFrameId = mCurrentFrameId;
                }
                
                projection += hnode->mCache.mProjection;
        }

        return projection;

#else

        if(node->mCache.mLastProjectionFrameId < mFrameId)
        {
                node->mCache.mProjection = 
                        ComputeScreenSpaceProjection(node->GetBox());
                node->mCache.mLastProjectionFrameId = mFrameId;
        }

        return node->mCache.mProjection;
#endif
}


float Bvh::ComputeScreenSpaceProjection(const BoundingBox &box) const
{
#if 0
        const float dist = CalcDistance(box);
#else
        //const float dist = sqrt(GetMinSquareDistance(box));
        const float dist = Distance(mCamera->GetPosition(), box.getCenter());
#endif

        const float scale = dist ? dist * mScale : 1e-6;
        const float f = 1.0f / scale;
        const float coverage = f * f * box.getSurface() / 6.0f;

#if 0
        OUT1("scale: " << scale);
        OUT1("box: " << box.getSurface());
        OUT1("aspect: " << aspect);
        OUT1("fov: " << fov);
        OUT1("tan: " << tan(fov * mypi / 180.0f));
        //OUT1("w: " << w << " h: " << h);
        OUT1("cov: " << coverage << " pix: " << coverage * w * h << " f: " << f);
#endif
        return coverage;
}


void Bvh::RenderBoundingBox(BvhNode *node)
{
#if 1
        static BvhNodeContainer dummy(1);
        dummy[0] = node;
        RenderBoundingBoxes(dummy);
#else
        RenderBoxIndividual(node);
#endif
}


void Bvh::RenderBoundingBoxImmediate(const BoundingBox &box, const bool restartStrip)
{
        const Point3f *pU = box.getUpperPtr();
        const Point3f *pL = box.getLowerPtr();

        ///////////
        //-- render AABB as triangle strips

        glVertex3f(pL->x, pL->y, pU->z);
        glVertex3f(pU->x, pL->y, pU->z);
        glVertex3f(pL->x, pU->y, pU->z);
        glVertex3f(pU->x, pU->y, pU->z);
        glVertex3f(pL->x, pU->y, pL->z); 
        glVertex3f(pU->x, pU->y, pL->z);
        glVertex3f(pL->x, pL->y, pL->z); 
        glVertex3f(pU->x, pL->y, pL->z);

        glPrimitiveRestartNV();

        //-- render second half of AABB
        glVertex3f(pL->x, pU->y, pU->z); 
        glVertex3f(pL->x, pU->y, pL->z);
        glVertex3f(pL->x, pL->y, pU->z);
        glVertex3f(pL->x, pL->y, pL->z);
        glVertex3f(pU->x, pL->y, pU->z); 
        glVertex3f(pU->x, pL->y, pL->z);
        glVertex3f(pU->x, pU->y, pU->z); 
        glVertex3f(pU->x, pU->y, pL->z);

        
        // restart for all but last node
        if (restartStrip)
        {
                glPrimitiveRestartNV();
        }
}


int Bvh::RenderBoundingBoxesImmediate(const BvhNodeContainer &nodes)
{
        timeSetup.Entry();
        int renderedBoxes = 0;
        
        glBegin(GL_TRIANGLE_STRIP);

        BvhNodeContainer::const_iterator bit, bit_end = nodes.end();

        for (bit = nodes.begin(); bit != bit_end; ++ bit)
        {
                BvhNode *node = *bit;

                renderedBoxes += node->mNumTestNodes;
#if 0
                const bool restartStrip = false;
                RenderBoundingBoxImmediate(node->GetBox(), restartStrip);
#else
                //OUT1("nodes: " << node->mTestNodes.size() << " l: " << node->IsLeaf());
                for (int i = 0; i < node->mNumTestNodes; ++ i)
                {
                        RenderableHierarchyNode *testNode = mTestNodes[node->mTestNodesIdx + i];
                        const BoundingBox &box = testNode->GetBox();

                        // restart for all but last node
                        const bool restartStrip = (bit != bit_end - 1) || (i != node->mNumTestNodes - 1);

                        RenderBoundingBoxImmediate(box, restartStrip);
                }
#endif
        }
        glEnd();
        timeSetup.Exit();

        return renderedBoxes;
}


int Bvh::RenderBoundingBoxes(const BvhNodeContainer &nodes)
{
        // always render in immediate mode if there is only one box to render
        if (((nodes.size() == 1) && (nodes[0]->mNumTestNodes == 1)) || 
                !Settings::Global()->get_nvocc_use_index_arrays())
        {       
                // render bounding boxes individually
                return RenderBoundingBoxesImmediate(nodes);
        }
        else 
        {
                const int renderedBoxes = 
                        PrepareBoundingBoxesWithDrawArrays(nodes);
                RenderBoundingBoxesWithDrawArrays(renderedBoxes);

                return renderedBoxes;
        }
}


int Bvh::PrepareBoundingBoxesWithDrawArrays(const BvhNodeContainer &nodes)
{
        //////
        //-- for the first time we come here ...

        if (!mIndices)
        {       // create list of indices
                CreateIndices();
        }

        if (sCurrentVboId == -1)
        {       
                // prepare the vbo 
                PrepareVertices();
        }

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

        timeSetup.Entry();

        int numNodes = 0;

        BvhNodeContainer::const_iterator nit, nit_end = nodes.end();

        for (nit = nodes.begin(); nit != nit_end; ++ nit)
        {
                BvhNode *node = *nit;

        const int numIndices = node->mNumTestNodes * sNumIndicesPerBox;
                
                // copy indices
                memcpy(mIndices + numNodes * sNumIndicesPerBox, 
                           mTestIndices + node->mIndicesPtr, 
#if 0 //ALIGN_INDICES
                           ((numIndices / 32) * 32 + 32) * sizeof(unsigned int));
#else
                           numIndices * sizeof(unsigned int));
#endif
                numNodes += node->mNumTestNodes;
        }

        timeSetup.Exit();

        return numNodes;
}


void Bvh::RenderBoundingBoxesWithDrawArrays(const int numNodes)
{
        //////
        //-- Rendering the vbo

        timeIssueDrawElements.Entry();

        if (useVbos)
                // set the vertex pointer to the vertex buffer
                glVertexPointer(3, GL_FLOAT, 0, (char *)NULL);          
        else
                glVertexPointer(3, GL_FLOAT, 0, mVertices);

        // we do use the last or the first index (they are generate and only used to connect strips)
        const int numElements = numNodes * sNumIndicesPerBox - 1;

        // don't render first degenerate index
        glDrawElements(GL_TRIANGLE_STRIP, numElements, GL_UNSIGNED_INT, mIndices + 1); 
        
        timeIssueDrawElements.Exit();
}


#define ALIGN_INDICES 1

void Bvh::CreateIndices()
{
        // collect bvh nodes
        BvhNodeContainer nodes;
        CollectNodes(mRoot, nodes);

        OUT1("creating new indices");

        int numMaxIndices = 0;

        BvhNodeContainer::const_iterator lit, lit_end = nodes.end();

        for (lit = nodes.begin(); lit != lit_end; ++ lit)
        {
                int offset = (*lit)->mNumTestNodes * sNumIndicesPerBox;
#if ALIGN_INDICES
                // align with 32
                offset = (offset / 32) * 32 + 32; 
#endif
                numMaxIndices += offset;
        }


        OUT1("creating global indices buffer");

        if (mIndices) delete [] mIndices;
        if (mTestIndices) delete [] mTestIndices;

        // global buffer: create it once so we don't have
        // to allocate memory from the chunks of the node
        mIndices = new unsigned int[numMaxIndices];

        // create new index buffer for the individual nodes
        mTestIndices = new unsigned int[numMaxIndices];
        
        mCurrentIndicesPtr = 0;

        for (lit = nodes.begin(); lit != lit_end; ++ lit)
        {
                BvhNode *node = *lit;
                
                // resize array
                node->mIndicesPtr = mCurrentIndicesPtr;

                int numIndices = 0;

                // the bounding box of the test nodes are rendered instead of the root node
                // => store their indices
                for (int i = 0; i < node->mNumTestNodes; ++ i, numIndices += sNumIndicesPerBox)
                {
                        RenderableHierarchyNode *testNode = mTestNodes[node->mTestNodesIdx + i];

                        // add indices to root node
                        for (int j = 0; j < sNumIndicesPerBox; ++ j)
                        {
                                mTestIndices[mCurrentIndicesPtr + numIndices + j] = sIndices[j] + testNode->GetId() * 8;
                        }
                }

                // align with 32
#if ALIGN_INDICES
                const int offset = (numIndices / 32) * 32 + 32;
#else
                const int offset = numIndices;
#endif
                mCurrentIndicesPtr += offset;
        }
}


void Bvh::ComputeIds()
{
        // collect all nodes, also the nodes from local bvh
        // warning: root nodes local bvh are not in there, as they 
        // are equivalent geometry bvh leaves
        HierarchyNodeContainer nodes;
        CollectAllNodes(mRoot, nodes);

        // assign ids to all nodes of the "regular" hierarchy
        int i = 0;
        HierarchyNodeContainer::const_iterator lit, lit_end = nodes.end();

        for (lit = nodes.begin(); lit != lit_end; ++ lit, ++ i)
        {
                (*lit)->SetId(i);
        }
}


void Bvh::PrepareVertices()
{
        // collect all nodes, also the nodes from local bvh
        // warning: root nodes local bvh are not in there, as they 
        // are equivalent geometry bvh leaves
        HierarchyNodeContainer nodes;

        nodes.reserve(GetNumNodes());
        CollectAllNodes(mRoot, nodes);

        const unsigned int bufferSize = 8 * (int)nodes.size();
        mVertices = new Point3f[bufferSize];
        
        int i = 0;
        
        // store bounding box vertices
        HierarchyNodeContainer::const_iterator lit, lit_end = nodes.end();

        for (lit = nodes.begin(); lit != lit_end; ++ lit, i += 8)
        {
                RenderableHierarchyNode *node = *lit;
                const Point3f *vertices = (const Point3f *)node->GetBox().getVertexData();

                for (int j = 0; j < 8; ++ j)
                {       
                        ((Point3f *)mVertices)[node->GetId() * 8 + j] = vertices[j];
                }
        }

        if (useVbos)
        {
                glGenBuffersARB(1, &sCurrentVboId);
                glBindBufferARB(GL_ARRAY_BUFFER_ARB, sCurrentVboId);
                glVertexPointer(3, GL_FLOAT, 0, (char *)NULL);
        
                glBufferDataARB(GL_ARRAY_BUFFER_ARB, 
                                    bufferSize * sizeof(Point3f), 
                                    mVertices, 
                                                GL_STATIC_DRAW_ARB);

                //glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);

                // data handled by graphics driver from now on
                DELAPTR(mVertices);

                OUT1("***** created vbos *********");
        }
        else
        {
                sCurrentVboId = 0;
        glVertexPointer(3, GL_FLOAT, 0, mVertices);

                OUT1("created vertices");
        }
}


void Bvh::SetMaxDepthForTestingChildren(const int maxDepth) 
{
        if (maxDepth != mMaxDepthForTestingChildren)
        {
                mMaxDepthForTestingChildren = maxDepth;
                RecomputeBounds();
        }
}


void Bvh::SetAreaRatioThresholdForTestingChildren(const float ratio) 
{
        if (ratio != mAreaRatioThreshold)
        {
                mAreaRatioThreshold = ratio;
                RecomputeBounds();
        }
}


void Bvh::SetVolRatioThresholdForTestingChildren(const float ratio) 
{
        if (ratio != mVolRatioThreshold)
        {
                mVolRatioThreshold = ratio;
                RecomputeBounds();
        }
}


void Bvh::SetUseTighterBoundsForTests(const bool tighterBoundsForTests)
{
        if (mUseTighterBoundsOnlyForLeafTests != tighterBoundsForTests)
        {
                mUseTighterBoundsOnlyForLeafTests = tighterBoundsForTests;
                RecomputeBounds();
        }
}


void Bvh::SetCollectTighterBoundsWithMaxLevel(bool t)
{
        if (mCollectTighterBoundsWithMaxLevel != t)
        {
                mCollectTighterBoundsWithMaxLevel = t;
                RecomputeBounds();
        }
}


void Bvh::RecomputeBounds()
{
        // clear old list
        mTestNodes.clear();

        // collect all nodes
        BvhNodeContainer nodes;
        CollectNodes(mRoot, nodes);

        OUT1("recomputing bounds, children will be tested in depth " << mMaxDepthForTestingChildren);

        if (mCollectTighterBoundsWithMaxLevel)
                OUT1("creating tighter bounds using max level");
        else
                OUT1("creating tighter bounds using best set");

        int success = 0;
        BvhNodeContainer::const_iterator lit, lit_end = nodes.end();

        for (lit = nodes.begin(); lit != lit_end; ++ lit)
        {
                BvhNode *node = *lit;

                if (mCollectTighterBoundsWithMaxLevel)
                {
                        //OUT1("creating tighter bounds using max level");

                        // recreate list of nodes that will be tested as a proxy ...
                        if (CreateNodeRenderList(node))
                                ++ success;
                }
                else
                {
                        //OUT1("creating tighter bounds using best set");

                        HierarchyNodeContainer hnodes;
                        CollectBestNodes(node, hnodes); 

                        // the new test nodes are added at the end of the vector
                        node->mTestNodesIdx = mTestNodes.size();

                        HierarchyNodeContainer::const_iterator cit;

                        // use the found nodes as nodes during the occlusion tests
                        for (cit = hnodes.begin(); cit != hnodes.end(); ++ cit)
                        {
                                RenderableHierarchyNode *child = *cit;
                                mTestNodes.push_back(child);
                        }

                        node->mNumTestNodes = (int)hnodes.size();
                }
                //OUT1("testnodes: " << node->mNumTestNodes);
        }

        const float p = 100.0f * (float)success / nodes.size();

        OUT1("created tighter bounds for " << p << " percent of the nodes");

        // recreate indices used for indirect mode rendering
        if (mIndices)
        {
                CreateIndices();
        }
}

        
bool Bvh::CreateNodeRenderList(BvhNode *node)
{
        HierarchyNodeContainer children;

        if (mUseTighterBoundsOnlyForLeafTests && !node->IsLeaf())
        {
                children.push_back(node);
        }
        else
        {
                // collect nodes that will be tested instead of the leaf node
                // in order to get a tighter bounding box test
                CollectNodes(node, mMaxDepthForTestingChildren, children);
        }


        // using the tighter bounds is not feasable in case
        // that the tighter bounds represent nearly the same projected area
        // as the old bounding box. Find this out using either volume or area
        // heuristics

        float vol = 0;
        float area = 0;

        HierarchyNodeContainer::const_iterator cit;

        for (cit = children.begin(); cit != children.end(); ++ cit)
        {
                RenderableHierarchyNode *c = *cit;

                area += c->GetBox().getSurface();
                vol += c->GetBox().getVolume();
        }

        const float volRatio = vol / node->GetBox().getVolume();
        const float areaRatio = area / node->GetBox().getSurface();
        
        bool success;

        if ((areaRatio < mAreaRatioThreshold) &&
                (volRatio < mVolRatioThreshold))
        {
                success = true;
        }
        else
        {
                // hack: only store node itself
                children.clear();
                children.push_back(node);

                success = false;
        }

        // the new test nodes are added at the end of the vector
        node->mTestNodesIdx = mTestNodes.size();

        // use the found nodes as nodes during the occlusion tests
        for (cit = children.begin(); cit != children.end(); ++ cit)
        {
                RenderableHierarchyNode *child = *cit;
                mTestNodes.push_back(child);
        }

        node->mNumTestNodes = (int)children.size();

        return success;
}


void Bvh::ResetNodeClassifications()
{
        BvhNodeContainer nodes;

        nodes.reserve(GetNumNodes());
        CollectNodes(mRoot, nodes);

        BvhNodeContainer::const_iterator lit, lit_end = nodes.end();

        for (lit = nodes.begin(); lit != lit_end; ++ lit)
        {
                (*lit)->ResetVisibility();
        }
}


int Bvh::PostProcessLeaves(BvhLeafContainer &leaves)
{
        OUT1("creating tighter bounds for leaves");

        BvhLeafContainer::const_iterator lit, lit_end = leaves.end();

        int i = 0;
        int tighter = 0;

        for (lit = leaves.begin(); lit != lit_end; ++ lit, ++ i)
        {
                BvhLeaf *leaf = *lit;

                int triangleCount;
                // collect the triangles from the stored geometry
                Point3f *triangles = CollectTriangles(leaf, triangleCount);

                if (CreateTighterBoundsForLeaf(leaf, triangles, triangleCount))
                        ++ tighter;

                leaf->mArea = ComputeGeometryArea(leaf, triangles, triangleCount);

                delete [] triangles;

                if (i % 1000 == 999)
                        OUT1(i << " leaves processed ");
        }

        const float p = 100.0f * tighter / leaves.size();
        OUT1("created tighter bounds for " << p << " percent of the leaves");

        return tighter;
}


bool Bvh::CreateTighterBoundsForLeaf(BvhLeaf *leaf, Point3f *triangles, const int triangleCount)
{
        // create a local bvh over the triangles
        if (leaf->mTriangleBvh) delete leaf->mTriangleBvh;

        leaf->mTriangleBvh = new TriangleBvh((Triangle *)triangles, triangleCount);

        const int maxDepth = Settings::Global()->get_nvocc_local_bvh_max_depth();
        const int maxTriangles = Settings::Global()->get_nvocc_local_bvh_max_triangles();
        const int minArea = Settings::Global()->get_nvocc_local_bvh_min_area();
        const int splitType = Settings::Global()->get_nvocc_local_bvh_split_type();

        leaf->mTriangleBvh->SetMaxDepth(maxDepth);
        leaf->mTriangleBvh->SetMaxTriangles(maxTriangles);
        leaf->mTriangleBvh->SetMinAreaRatio(minArea);
        leaf->mTriangleBvh->SetSplitType(splitType);
        
        leaf->mTriangleBvh->SetVerboseOutput(false);

        // construct
        leaf->mTriangleBvh->Construct();

        return true;
}


float Bvh::ComputeGeometryArea(BvhLeaf *leaf, Point3f *triangles, const int triangleCount) const
{
        float area = 0;

        for (int i = 0; i < triangleCount; ++ i)
        {
                area += ((Triangle *)triangles)[i].GetArea();
        }

        return area;
}


Point3f *Bvh::CollectTriangles(BvhLeaf *leaf, int &triangleCount)
{
        // count triangles in the leaf
        triangleCount = 0;

        for (int geomIdx = leaf->mFirst; geomIdx <= leaf->mLast; ++ geomIdx)
        {
                Shape3D *sh = static_cast<Shape3D *>(mGeometry[geomIdx]->GetNode());
                IndexedTriangleStripArray *strips = static_cast<IndexedTriangleStripArray *>(sh->getGeometry());

                triangleCount += strips->getNumTriangles();
        }

        Triangle *triangles = new Triangle[triangleCount];
        //OUT1("The leaf contains " << triangleCount << " triangles in " << (int)geometry.size() << " nodes");

        int currentIdx = 0;
        
        for (int geomIdx = leaf->mFirst; geomIdx <= leaf->mLast; ++ geomIdx)
        {
                NodeGeometry *geom = mGeometry[geomIdx];

                Shape3D *sh = static_cast<Shape3D *>(geom->GetNode());
                
                IndexedTriangleStripArray *strips = 
                        static_cast<IndexedTriangleStripArray *>(sh->getGeometry());

                // get indices of containted triangles (strips are converted before)
                int *indices;
        
                int tCount = strips->getNumTriangles();
                int indexCount = tCount * 3;
                
                // compute triangles from strips
                Point3f *coordinates = (Point3f *)strips->getCoordRef3f();
                strips->getTriangleCoordinateIndices(0, &indices, indexCount);

                for (int i = 0; i < tCount; ++ i)
                {
                        triangles[i + currentIdx].mVertices[0] = coordinates[indices[i * 3 + 0]]; 
                        triangles[i + currentIdx].mVertices[1] = coordinates[indices[i * 3 + 1]]; 
                        triangles[i + currentIdx].mVertices[2] = coordinates[indices[i * 3 + 2]]; 

                        if (geom->mMatId != NV_RenderPredictorRenderAction::NO_MAT)
                        {
                                const Matrix4f &m = mRenderer->app_mats[geom->mMatId];
                                Transform3D tf(m);
                                tf.transform(triangles[i + currentIdx].mVertices[0]);
                                tf.transform(triangles[i + currentIdx].mVertices[1]);
                                tf.transform(triangles[i + currentIdx].mVertices[2]);
                        }
                }

                currentIdx += tCount;
        }

        return (Point3f *)triangles;
}


int Bvh::CountTriangles(BvhLeaf *leaf) const
{
        int triangleCount = 0;
        
        for (int i = leaf->mFirst; i <= leaf->mLast; ++ i)
        {
                triangleCount += mGeometry[i]->mNumTriangles;
        }

        return triangleCount;
}


void Bvh::ComputeBvhStats() 
{
        std::stack<BvhNode *> nodeStack;
        nodeStack.push(mRoot);

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

                if (node->IsLeaf()) 
                {
                        BvhLeaf *leaf = static_cast<BvhLeaf *>(node);

                        mBvhStats.mLeafSA += leaf->mBox.getSurface();
                        mBvhStats.mLeafVol += leaf->mBox.getVolume();

                        TriangleBvh::TriangleBvhStats tStats;
                        leaf->mTriangleBvh->GetBvhStats(tStats);

                        mBvhStats.mBoundsLeafSA += tStats.mLeafSA;
                        mBvhStats.mBoundsInteriorSA += tStats.mInteriorSA;

                        mBvhStats.mBoundsLeafVol += tStats.mLeafVol;
                        mBvhStats.mBoundsInteriorVol += tStats.mInteriorVol;

                        mBvhStats.mBoundsLeavesCount += leaf->mTriangleBvh->GetNumLeaves();
                } 
                else 
                {
                        mBvhStats.mInteriorSA += node->mBox.getSurface();
                        mBvhStats.mInteriorVol += node->mBox.getVolume();

                        BvhInterior *interior = static_cast<BvhInterior *>(node);
                        
                        nodeStack.push(interior->mBack);
                        nodeStack.push(interior->mFront);
                }
        }

        mBvhStats.mGeometryRatio = mGeometrySize / (float)GetNumLeaves();
        mBvhStats.mTriangleRatio = mBvhStats.mTriangles / (float)GetNumLeaves();
}


void Bvh::PrintBvhStats() const
{
        //OUT1("triangle bvh surface " << mBvhStats.mBoundsLeafSA);
        //OUT1("triangle bvh volume " << mBvhStats.mBoundsLeafVol);

        OUT1("bvh stats:");
        OUT1("interiorNodesSA = " << mBvhStats.mInteriorSA / mRoot->mBox.getSurface());
        OUT1("leafNodesSA = " << mBvhStats.mLeafSA / mRoot->mBox.getSurface());
        OUT1("interiorNodesVolume = " << mBvhStats.mInteriorVol / mRoot->mBox.getVolume());
        OUT1("leafNodesVolume = " << mBvhStats.mLeafVol / mRoot->mBox.getVolume() << "\n");

        OUT1("boundsInteriorNodesSA = " << mBvhStats.mBoundsInteriorSA / mBvhStats.mLeafSA);
        OUT1("boundsLeafNodesSA = " << mBvhStats.mBoundsLeafSA / mBvhStats.mLeafSA);
        OUT1("boundsInteriorNodesVolume = " << mBvhStats.mBoundsInteriorVol / mBvhStats.mLeafVol);
        OUT1("boundsLeafNodesVolume = " << mBvhStats.mBoundsLeafVol / mBvhStats.mLeafVol << "\n");

        OUT1("boundsLeaves: " <<  (float)mBvhStats.mBoundsLeavesCount / (float)GetNumLeaves() << "\n");
        OUT1("geometry per leaf: " <<  mBvhStats.mGeometryRatio);
        OUT1("triangles per leaf: " <<  mBvhStats.mTriangleRatio);
        OUT1("");
}


static void RenderBoxForViz(const BoundingBox &box)
{
        glBegin(GL_LINE_LOOP);
        glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
        glEnd();

        glBegin(GL_LINE_LOOP);
        glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
        glEnd();

        glBegin(GL_LINE_LOOP);
        glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
        glEnd();

        glBegin(GL_LINE_LOOP);
        glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
        glEnd();

        glBegin(GL_LINE_LOOP);
        glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
        glEnd();

        glBegin(GL_LINE_LOOP);
        glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
        glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
        glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);

        glEnd();
}


void Bvh::RenderBoundingBoxesForViz(const int mode)
{
        BvhLeafContainer leaves;
        leaves.reserve(mRoot->GetNumLeaves());

        CollectLeaves(mRoot, leaves);

        BvhLeafContainer::const_iterator it, it_end = leaves.end();

        for (it = leaves.begin(); it != it_end; ++ it)
        {
                BvhLeaf *leaf = *it;

                if ((mode == 0) || (mode == 2))
                {
                        glColor3f(1.0f, 1.0f, 1.0f);
            RenderBoxForViz(leaf->GetBox());
                }
                
                if ((mode == 1) || (mode == 2)) 
                {
                        glColor3f(1.0f, 0, 0);
                        for (int i = 0; i < leaf->mNumTestNodes; ++ i)
                        {
                                RenderBoxForViz(mTestNodes[leaf->mTestNodesIdx + i]->GetBox());
                        }
                }
        }
}


int Bvh::CountTriangles(BvhNode *node) const
{
        int numTriangles = 0;

        for (int i = node->mFirst; i <= node->mLast; ++ i)
        {
                numTriangles += mGeometry[i]->mNumTriangles;
        }

        return numTriangles;
}


float Bvh::GetGeometryArea(BvhNode *node) const
{
        return node->mArea;
}


float Bvh::GetAvgDepth() const
{
        return mAvgDepth;
}


static float GetNodePriority(RenderableHierarchyNode *node)
{
        if (node->IsLeaf())
                return 0.0f;

        float result;

        if (node->GetType() == BVH_NODE)
        {
                BvhInterior *interior = static_cast<BvhInterior *>(node);
                
                // evaluate the priority of this node
                const float correctionFactor = 0.0f;
                //  1.0f*interior->box.FaceArea(node->axis);

                result = 
                        interior->GetBox().getSurface() - 
                        (interior->GetBack()->GetBox().getSurface() + interior->GetFront()->GetBox().getSurface()) + 
                        correctionFactor;
        }
        else
        {
                TriangleBvhInterior *interior = static_cast<TriangleBvhInterior *>(node);
                
                // evaluate the priority of this node
                const float correctionFactor = 0.0f;
                //  1.0f*interior->box.FaceArea(node->axis);

                result = 
                        interior->GetBox().getSurface() - 
                        (interior->GetBack()->GetBox().getSurface() + interior->GetFront()->GetBox().getSurface()) + 
                        correctionFactor;
        }

        return result;
}


struct NodeEntry
{
        NodeEntry(RenderableHierarchyNode *node, float p): mNode(node), mPriority(p)
        {}

        RenderableHierarchyNode *mNode;
        float mPriority;
};


bool operator<(const NodeEntry &a, const NodeEntry &b)
{
        return a.mPriority < b.mPriority;
}


void Bvh::CollectBestNodes(RenderableHierarchyNode *node, HierarchyNodeContainer &nodes)
{
        int maxNodes = 32;
        priority_queue<NodeEntry> nodeStack;

        nodeStack.push(NodeEntry(node, GetNodePriority(node)));
        nodes.clear();

        float SA = node->GetBox().getSurface();
        float saThreshold = 1.1f * SA;
        int numNodes = 1;
        OUT1(numNodes << " " << SA);

        while (!nodeStack.empty() &&
                int(nodes.size() + nodeStack.size()) < maxNodes)
        {
                NodeEntry entry = nodeStack.top();
                nodeStack.pop();

                RenderableHierarchyNode *current = entry.mNode;

                if (current->IsLeaf())
                {
                        // check if there further subdivision on triangle level
                        if ((current->GetType() == BVH_NODE) && (((BvhLeaf *)current)->mTriangleBvh->GetNumNodes() > 1))
                        {
                                // push back the root of the local bvh
                                nodeStack.push(NodeEntry(((BvhLeaf *)current)->mTriangleBvh->GetRoot(), entry.mPriority));
                        }
                        else // terminate traversal
                        {
                                nodes.push_back(current);
                        }
                }
                else // interior node
                {
                        // surface area of child nodes
                        SA -= entry.mPriority;
                
                        // finish the search if we have too much fillrate increase
                        if (SA > saThreshold)
                        {
                                OUT1("break at " << SA << " > " << saThreshold);
                                nodes.push_back(current);
                        }
                        else
                        {
                                ++ numNodes;

                                OUT1(numNodes << " " << SA << " " << entry.mPriority);

                                if (current->GetType() == BVH_NODE)
                                {
                                        BvhInterior *interior = static_cast<BvhInterior *>(current);

                                        nodeStack.push(NodeEntry(interior->GetBack(), 
                                                GetNodePriority(interior->GetBack())));
                                        nodeStack.push(NodeEntry(interior->GetFront(),
                                                GetNodePriority(interior->GetFront())));
                                }
                                else
                                {
                                        TriangleBvhInterior *interior = static_cast<TriangleBvhInterior *>(current);

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

        while (!nodeStack.empty()) 
        {
                NodeEntry entry = nodeStack.top();
        
                nodeStack.pop();
                RenderableHierarchyNode *current = entry.mNode;
                nodes.push_back(current);
        }
} 

#endif