source: GTP/trunk/Lib/Vis/Preprocessing/src/BvHierarchy.cpp @ 2255

Revision 2255, 86.2 KB checked in by mattausch, 18 years ago (diff)

improved scenemanager config

Line 
1#include <stack>
2#include <time.h>
3#include <iomanip>
4
5#include "BvHierarchy.h"
6#include "ViewCell.h"
7#include "Plane3.h"
8#include "Mesh.h"
9#include "common.h"
10#include "Environment.h"
11#include "Polygon3.h"
12#include "Ray.h"
13#include "AxisAlignedBox3.h"
14#include "Exporter.h"
15#include "Plane3.h"
16#include "ViewCellsManager.h"
17#include "Beam.h"
18#include "VspTree.h"
19#include "HierarchyManager.h"
20
21
22namespace GtpVisibilityPreprocessor {
23
24
25#define PROBABILIY_IS_BV_VOLUME 1
26#define USE_FIXEDPOINT_T 0
27#define USE_VOLUMES_FOR_HEURISTICS 1
28#define TEST_POWERPLANT 0
29#define USE_BETTER_RENDERCOST_EST 0
30
31//int BvhNode::sMailId = 10000;
32//int BvhNode::sReservedMailboxes = 1;
33
34BvHierarchy *BvHierarchy::BvhSubdivisionCandidate::sBvHierarchy = NULL;
35
36
37/// sorting operator
38inline static bool ilt(Intersectable *obj1, Intersectable *obj2)
39{
40        return obj1->mId < obj2->mId;
41}
42
43
44/// sorting operator
45inline static bool smallerSize(Intersectable *obj1, Intersectable *obj2)
46{
47        return obj1->GetBox().SurfaceArea() < obj2->GetBox().SurfaceArea();
48}
49
50
51
52/***************************************************************/
53/*              class BvhNode implementation                   */
54/***************************************************************/
55
56
57BvhNode::BvhNode():
58mParent(NULL),
59mTimeStamp(0),
60mRenderCost(-1)
61
62{
63       
64}
65
66BvhNode::BvhNode(const AxisAlignedBox3 &bbox):
67mParent(NULL),
68mBoundingBox(bbox),
69mTimeStamp(0),
70mRenderCost(-1)
71{
72}
73
74
75BvhNode::BvhNode(const AxisAlignedBox3 &bbox, BvhInterior *parent):
76mBoundingBox(bbox),
77mParent(parent),
78mTimeStamp(0),
79mRenderCost(-1)
80{
81}
82
83
84bool BvhNode::IsRoot() const
85{
86        return mParent == NULL;
87}
88
89
90BvhInterior *BvhNode::GetParent()
91{
92        return mParent;
93}
94
95
96void BvhNode::SetParent(BvhInterior *parent)
97{
98        mParent = parent;
99}
100
101
102int BvhNode::GetRandomEdgePoint(Vector3 &point,
103                                                                Vector3 &normal)
104{
105        // get random edge
106        const int idx = Random(12);
107        Vector3 a, b;
108        mBoundingBox.GetEdge(idx, &a, &b);
109       
110        const float w = RandomValue(0.0f, 1.0f);
111
112        point = a * w + b * (1.0f - w);
113
114        // TODO
115        normal = Vector3(0);
116
117        return idx;
118}
119
120
121
122/******************************************************************/
123/*              class BvhInterior implementation                  */
124/******************************************************************/
125
126
127BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox):
128BvhNode(bbox),
129mSubdivisionCandidate(NULL),
130mGlList(0)
131{
132  mActiveNode = this;
133}
134
135
136BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox, BvhInterior *parent):
137  BvhNode(bbox, parent),
138  mGlList(0)
139 
140{
141        mActiveNode = this;
142}
143
144
145BvhLeaf::BvhLeaf(const AxisAlignedBox3 &bbox,
146                                 BvhInterior *parent,
147                                 const int numObjects):
148  BvhNode(bbox, parent),
149  mGlList(0)
150
151{
152        mObjects.reserve(numObjects);
153        mActiveNode = this;
154}
155
156
157bool BvhLeaf::IsLeaf() const
158{
159        return true;
160}
161
162
163BvhLeaf::~BvhLeaf()
164{
165}
166
167
168void BvhLeaf::CollectObjects(ObjectContainer &objects)
169{
170        ObjectContainer::const_iterator oit, oit_end = mObjects.end();
171        for (oit = mObjects.begin(); oit != oit_end; ++ oit)
172        {
173                objects.push_back(*oit);
174        }
175}
176
177
178
179/******************************************************************/
180/*              class BvhInterior implementation                  */
181/******************************************************************/
182
183
184BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox):
185BvhNode(bbox), mFront(NULL), mBack(NULL)
186{
187}
188
189
190BvhInterior::BvhInterior(const AxisAlignedBox3 &bbox, BvhInterior *parent):
191BvhNode(bbox, parent), mFront(NULL), mBack(NULL)
192{
193}
194
195
196void BvhInterior::ReplaceChildLink(BvhNode *oldChild, BvhNode *newChild)
197{
198        if (mBack == oldChild)
199                mBack = newChild;
200        else
201                mFront = newChild;
202}
203
204
205bool BvhInterior::IsLeaf() const
206{
207        return false;
208}
209
210
211BvhInterior::~BvhInterior()
212{
213        DEL_PTR(mFront);
214        DEL_PTR(mBack);
215}
216
217
218void BvhInterior::SetupChildLinks(BvhNode *front, BvhNode *back)
219{
220    mBack = back;
221    mFront = front;
222}
223
224
225void BvhInterior::CollectObjects(ObjectContainer &objects)
226{
227        mFront->CollectObjects(objects);
228        mBack->CollectObjects(objects);
229}
230
231
232/*******************************************************************/
233/*                  class BvHierarchy implementation               */
234/*******************************************************************/
235
236
237BvHierarchy::BvHierarchy():
238mRoot(NULL),
239mTimeStamp(1),
240mIsInitialSubdivision(false)
241{
242        ReadEnvironment();
243        mSubdivisionCandidates = new SortableEntryContainer;
244//      for (int i = 0; i < 4; ++ i)
245//              mSortedObjects[i] = NULL;
246}
247
248
249BvHierarchy::~BvHierarchy()
250{
251        // delete the local subdivision candidates
252        DEL_PTR(mSubdivisionCandidates);
253
254        // delete the presorted objects
255        /*for (int i = 0; i < 4; ++ i)
256        {
257                DEL_PTR(mSortedObjects[i]);
258        }*/
259       
260        // delete the tree
261        DEL_PTR(mRoot);
262}
263
264
265void BvHierarchy::ReadEnvironment()
266{
267        bool randomize = false;
268        Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.randomize", randomize);
269         
270        // initialise random generator for heuristics
271        if (randomize)
272                Randomize();
273
274        //////////////////////////////
275        //-- termination criteria for autopartition
276
277        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.maxDepth", mTermMaxDepth);
278        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.maxLeaves", mTermMaxLeaves);
279        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.minObjects", mTermMinObjects);
280        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.minRays", mTermMinRays);
281        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.minProbability", mTermMinProbability);
282    Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.missTolerance", mTermMissTolerance);
283
284
285        //////////////////////////////
286        //-- max cost ratio for early tree termination
287
288        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.maxCostRatio", mTermMaxCostRatio);
289        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Termination.minGlobalCostRatio",
290                mTermMinGlobalCostRatio);
291        Environment::GetSingleton()->GetIntValue("BvHierarchy.Termination.globalCostMissTolerance",
292                mTermGlobalCostMissTolerance);
293
294
295        //////////////////////////////
296        //-- factors for subdivision heuristics
297
298        // if only the driving axis is used for splits
299        Environment::GetSingleton()->GetBoolValue("BvHierarchy.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
300        Environment::GetSingleton()->GetFloatValue("BvHierarchy.maxStaticMemory", mMaxMemory);
301        Environment::GetSingleton()->GetBoolValue("BvHierarchy.useCostHeuristics", mUseCostHeuristics);
302        Environment::GetSingleton()->GetBoolValue("BvHierarchy.useSah", mUseSah);
303
304    char subdivisionStatsLog[100];
305        Environment::GetSingleton()->GetStringValue("BvHierarchy.subdivisionStats", subdivisionStatsLog);
306        mSubdivisionStats.open(subdivisionStatsLog);
307
308        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
309        Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.useGlobalSorting", mUseGlobalSorting);
310        Environment::GetSingleton()->GetIntValue("BvHierarchy.minRaysForVisibility", mMinRaysForVisibility);
311        Environment::GetSingleton()->GetIntValue("BvHierarchy.maxTests", mMaxTests);
312        Environment::GetSingleton()->GetBoolValue("BvHierarchy.Construction.useInitialSubdivision", mApplyInitialPartition);
313        Environment::GetSingleton()->GetIntValue("BvHierarchy.Construction.Initial.minObjects", mInitialMinObjects);
314        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.Initial.maxAreaRatio", mInitialMaxAreaRatio);
315        Environment::GetSingleton()->GetFloatValue("BvHierarchy.Construction.Initial.minArea", mInitialMinArea);
316
317        //mMemoryConst = (float)(sizeof(VspLeaf) + sizeof(VspViewCell));
318        //mMemoryConst = (float)sizeof(BvhLeaf);
319        mMemoryConst = 16;//(float)sizeof(ObjectContainer);
320
321    mUseBboxAreaForSah = true;
322
323        /////////////
324        //-- debug output
325
326        Debug << "******* Bvh hierarchy options ******** " << endl;
327    Debug << "max depth: " << mTermMaxDepth << endl;
328        Debug << "min probabiliy: " << mTermMinProbability<< endl;
329        Debug << "min objects: " << mTermMinObjects << endl;
330        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
331        Debug << "miss tolerance: " << mTermMissTolerance << endl;
332        Debug << "max leaves: " << mTermMaxLeaves << endl;
333        Debug << "randomize: " << randomize << endl;
334        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
335        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
336        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
337        Debug << "max memory: " << mMaxMemory << endl;
338        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
339        Debug << "use surface area heuristics: " << mUseSah << endl;
340        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
341        //Debug << "split borders: " << mSplitBorder << endl;
342        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
343        Debug << "use global sort: " << mUseGlobalSorting << endl;
344        Debug << "minimal rays for visibility: " << mMinRaysForVisibility << endl;
345        Debug << "bvh mem const: " << mMemoryConst << endl;
346        Debug << "apply initial partition: " << mApplyInitialPartition << endl;
347        Debug << "min objects: " << mInitialMinObjects << endl;
348        Debug << "max area ratio: " << mInitialMaxAreaRatio << endl;
349        Debug << "min area: " << mInitialMinArea << endl;
350
351        Debug << endl;
352}
353
354
355void BvHierarchy::AssociateObjectsWithLeaf(BvhLeaf *leaf)
356{
357        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
358
359        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
360        {
361                (*oit)->mBvhLeaf = leaf;
362        }
363}
364
365
366static int CountRays(const ObjectContainer &objects)
367{
368        int nRays = 0;
369
370        ObjectContainer::const_iterator oit, oit_end = objects.end();
371
372        // warning: not exact number (there can be rays counted twice)
373        // otherwise we would have to traverse trough all rays
374        for (oit = objects.begin(); oit != oit_end; ++ oit)
375        {
376                nRays += (int)(*oit)->GetOrCreateRays()->size();
377        }
378
379        return nRays;
380}
381
382                                                                       
383float BvHierarchy::GetViewSpaceVolume() const
384{
385        return mViewCellsManager->GetViewSpaceBox().GetVolume();
386}
387
388
389BvhInterior *BvHierarchy::SubdivideNode(const BvhSubdivisionCandidate &sc,
390                                                                                BvhTraversalData &frontData,
391                                                                                BvhTraversalData &backData)
392{
393        // fill view cells cache
394        mNodeTimer.Entry();
395
396        const BvhTraversalData &tData = sc.mParentData;
397        BvhLeaf *leaf = tData.mNode;
398
399        AxisAlignedBox3 parentBox = leaf->GetBoundingBox();
400
401        // update stats: we have two new leaves
402        mBvhStats.nodes += 2;
403
404        if (tData.mDepth > mBvhStats.maxDepth)
405        {
406                mBvhStats.maxDepth = tData.mDepth;
407        }
408
409        // add the new nodes to the tree
410        BvhInterior *node = new BvhInterior(parentBox, leaf->GetParent());
411
412
413        //////////////////
414        //-- create front and back leaf
415
416        AxisAlignedBox3 fbox = EvalBoundingBox(sc.mFrontObjects, &parentBox);
417        AxisAlignedBox3 bbox = EvalBoundingBox(sc.mBackObjects, &parentBox);
418
419        BvhLeaf *back = new BvhLeaf(bbox, node, (int)sc.mBackObjects.size());
420        BvhLeaf *front = new BvhLeaf(fbox, node, (int)sc.mFrontObjects.size());
421
422        BvhInterior *parent = leaf->GetParent();
423
424        // replace a link from node's parent
425        if (parent)
426        {
427                parent->ReplaceChildLink(leaf, node);
428                node->SetParent(parent);
429        }
430        else // no parent => this node is the root
431        {
432                mRoot = node;
433        }
434
435        // and setup child links
436        node->SetupChildLinks(front, back);
437
438        ++ mBvhStats.splits;
439
440 
441        ////////////////////////////////////////
442        //-- fill front and back traversal data with the new values
443
444        frontData.mDepth = backData.mDepth = tData.mDepth + 1;
445
446        frontData.mNode = front;
447        backData.mNode = back;
448
449        backData.mSampledObjects = new ObjectContainer();
450        frontData.mSampledObjects = new ObjectContainer();
451
452        *backData.mSampledObjects = sc.mSampledBackObjects;
453        *frontData.mSampledObjects = sc.mSampledFrontObjects;
454
455        back->mObjects = sc.mBackObjects;
456        front->mObjects = sc.mFrontObjects;
457
458        // if the number of rays is too low, no assumptions can be made
459        // (=> switch to surface area heuristics?)
460        frontData.mNumRays = CountRays(sc.mSampledFrontObjects);
461        backData.mNumRays = CountRays(sc.mSampledBackObjects);
462
463        AssociateObjectsWithLeaf(back);
464        AssociateObjectsWithLeaf(front);
465 
466        ////////////
467        //-- compute pvs correction to cope with undersampling
468
469        frontData.mPvs = (float)sc.mNumFrontViewCells;//(float)CountViewCells(sc.mSampledFrontObjects);
470        backData.mPvs = (float)sc.mNumBackViewCells;//(float)CountViewCells(sc.mSampledBackObjects);
471
472        frontData.mCorrectedPvs = sc.mCorrectedFrontPvs;
473        backData.mCorrectedPvs = sc.mCorrectedBackPvs;
474
475
476        // compute probability of this node being visible,
477        // i.e., volume of the view cells that can see this node
478        frontData.mVolume = sc.mVolumeFrontViewCells;//EvalViewCellsVolume(sc.mSampledFrontObjects) / GetViewSpaceVolume();
479        backData.mVolume = sc.mVolumeBackViewCells;//EvalViewCellsVolume(sc.mSampledBackObjects) / GetViewSpaceVolume();
480
481        frontData.mCorrectedVolume = sc.mCorrectedFrontVolume;
482        backData.mCorrectedVolume = sc.mCorrectedBackVolume;
483
484
485    // how often was max cost ratio missed in this branch?
486        frontData.mMaxCostMisses = sc.GetMaxCostMisses();
487        backData.mMaxCostMisses = sc.GetMaxCostMisses();
488
489        // set the time stamp so the order of traversal can be reconstructed
490        node->SetTimeStamp(mHierarchyManager->mTimeStamp ++);
491         
492        // assign the objects in sorted order
493        if (mUseGlobalSorting)
494        {
495                AssignSortedObjects(sc, frontData, backData);
496        }
497
498        mNodeTimer.Exit();
499
500        // return the new interior node
501        return node;
502}
503
504
505BvhNode *BvHierarchy::Subdivide(SplitQueue &tQueue,
506                                                                SubdivisionCandidate *splitCandidate,
507                                                                const bool globalCriteriaMet
508                                                                ,vector<SubdivisionCandidate *> &dirtyList
509                                                                )
510{
511        mSubdivTimer.Entry();
512
513        BvhSubdivisionCandidate *sc =
514                static_cast<BvhSubdivisionCandidate *>(splitCandidate);
515        BvhTraversalData &tData = sc->mParentData;
516
517        BvhNode *currentNode = tData.mNode;
518
519        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
520        {       
521                //////////////
522                //-- continue subdivision
523
524                BvhTraversalData tFrontData;
525                BvhTraversalData tBackData;
526               
527                // create new interior node and two leaf node
528                currentNode = SubdivideNode(*sc, tFrontData, tBackData);
529
530                // decrease the weighted average cost of the subdivisoin
531                mTotalCost -= sc->GetRenderCostDecrease();
532                mPvsEntries += sc->GetPvsEntriesIncr();
533
534                // subdivision statistics
535                if (1) PrintSubdivisionStats(*sc);
536
537
538                ///////////////////////////
539                //-- push the new split candidates on the queue
540               
541                BvhSubdivisionCandidate *frontCandidate =
542                                new BvhSubdivisionCandidate(tFrontData);
543                BvhSubdivisionCandidate *backCandidate =
544                                new BvhSubdivisionCandidate(tBackData);
545               
546                // preprocess view cells
547                AssociateViewCellsWithObjects(*tData.mSampledObjects);
548
549                EvalSubdivisionCandidate(*frontCandidate, true, false);
550                EvalSubdivisionCandidate(*backCandidate, true, false);
551
552                CollectDirtyCandidates(sc, dirtyList, true);
553               
554                // release preprocessed view cells
555                ReleaseViewCells(*tData.mSampledObjects);
556               
557                // cross reference
558                tFrontData.mNode->SetSubdivisionCandidate(frontCandidate);
559                tBackData.mNode->SetSubdivisionCandidate(backCandidate);
560
561                //cout << "f: " << frontCandidate->GetPriority() << " b: " << backCandidate->GetPriority() << endl;
562                tQueue.Push(frontCandidate);
563                tQueue.Push(backCandidate);
564        }
565
566        /////////////////////////////////
567        //-- node is a leaf => terminate traversal
568
569        if (currentNode->IsLeaf())
570        {
571                /////////////////////
572                //-- store additional info
573                EvaluateLeafStats(tData);
574       
575                // this leaf is no candidate for splitting anymore
576                // => detach subdivision candidate
577                tData.mNode->SetSubdivisionCandidate(NULL);
578                // detach node so we don't delete it with the traversal data
579                tData.mNode = NULL;
580        }
581       
582        mSubdivTimer.Exit();
583
584        return currentNode;
585}
586
587
588float BvHierarchy::EvalPriority(const BvhSubdivisionCandidate &splitCandidate,
589                                                                const float renderCostDecr,
590                                                                const float oldRenderCost) const
591{
592        float priority;
593
594        if (mIsInitialSubdivision)
595        {
596                priority = (float)-splitCandidate.mParentData.mDepth;
597                return priority;
598        }
599
600        BvhLeaf *leaf = splitCandidate.mParentData.mNode;
601
602        // surface area heuristics is used when there is
603        // no view space subdivision available.
604        // In order to have some prioritized traversal,
605        // we use this formula instead
606        if (mHierarchyManager->GetViewSpaceSubdivisionType() ==
607                HierarchyManager::NO_VIEWSPACE_SUBDIV)
608        {
609                priority = EvalSahCost(leaf);
610        }
611        else
612        {
613                // take render cost of node into account
614                // otherwise danger of being stuck in a local minimum!
615                const float factor = mRenderCostDecreaseWeight;
616
617                priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;
618               
619                if (mHierarchyManager->mConsiderMemory)
620                {
621                        priority /= ((float)splitCandidate.GetPvsEntriesIncr() + mMemoryConst);
622                }
623        }
624
625        // hack: don't allow empty splits to be taken
626        if (splitCandidate.mFrontObjects.empty() || splitCandidate.mBackObjects.empty())
627                priority = 0;
628
629        return priority;
630}
631
632
633static float AvgRayContribution(const int pvs, const int nRays)
634{
635        return (float)pvs / ((float)nRays + Limits::Small);
636}
637
638
639static float AvgRaysPerObject(const int pvs, const int nRays)
640{
641        return (float)nRays / ((float)pvs + Limits::Small);
642}
643
644
645void BvHierarchy::EvalSubdivisionCandidate(BvhSubdivisionCandidate &splitCandidate,
646                                                                                   const bool computeSplitPlane,
647                                                                                   const bool preprocessViewCells)
648{
649        mPlaneTimer.Entry();
650
651#if STORE_VIEWCELLS_WITH_BVH
652        if (preprocessViewCells) // fill view cells cache
653                AssociateViewCellsWithObjects(*splitCandidate.mParentData.mSampledObjects);
654#endif
655
656        if (computeSplitPlane)
657        {
658                splitCandidate.mFrontObjects.clear();
659                splitCandidate.mBackObjects.clear();
660                splitCandidate.mSampledFrontObjects.clear();
661                splitCandidate.mSampledBackObjects.clear();
662
663                const bool sufficientSamples =
664                        splitCandidate.mParentData.mNumRays > mMinRaysForVisibility;
665
666                //if (!sufficientSamples) cout << splitCandidate.mParentData.mNumRays << " ";
667
668                const bool useVisibiliyBasedHeuristics =
669                                        !mUseSah &&
670                                        (mHierarchyManager->GetViewSpaceSubdivisionType() ==
671                                        HierarchyManager::KD_BASED_VIEWSPACE_SUBDIV) &&
672                                        sufficientSamples;
673
674                // compute best object partition
675                const float ratio =     SelectObjectPartition(splitCandidate.mParentData,
676                                                                                                  splitCandidate.mFrontObjects,
677                                                                                                  splitCandidate.mBackObjects,
678                                                                                                  useVisibiliyBasedHeuristics);
679       
680                // cost ratio violated?
681                const bool maxCostRatioViolated = mTermMaxCostRatio < ratio;
682                const int previousMisses = splitCandidate.mParentData.mMaxCostMisses;
683
684                splitCandidate.SetMaxCostMisses(maxCostRatioViolated ?
685                                                                                previousMisses + 1 : previousMisses);
686
687                StoreSampledObjects(splitCandidate.mSampledFrontObjects, splitCandidate.mFrontObjects);
688                StoreSampledObjects(splitCandidate.mSampledBackObjects, splitCandidate.mBackObjects);
689        }
690
691        mPlaneTimer.Exit();
692
693        mEvalTimer.Entry();
694
695        const BvhTraversalData &tData = splitCandidate.mParentData;
696        BvhLeaf *leaf = tData.mNode;
697
698        // avg contribution of a ray to a pvs
699        const float pvs = (float)CountViewCells(*tData.mSampledObjects);
700        //const float avgRayContri = AvgRayContribution((int)pvs, tData.mNumRays);
701        const float avgRaysPerObject = AvgRaysPerObject((int)pvs, tData.mNumRays);
702
703        // high avg ray contri, the result is influenced by undersampling
704        splitCandidate.SetAvgRaysPerObject(avgRaysPerObject);
705        const float viewSpaceVol = GetViewSpaceVolume();
706
707        const float oldVolume = EvalViewCellsVolume(*tData.mSampledObjects) / viewSpaceVol;
708        const float oldRatio = (tData.mVolume > 0) ? oldVolume / tData.mVolume : 1;
709        const float parentVol = tData.mCorrectedVolume * oldRatio;
710
711        //cout << "h " << avgRaysPerObject << " ";
712        // this leaf is a pvs entry in all the view cells
713        // that see one of the objects.
714        splitCandidate.mVolumeFrontViewCells = EvalViewCellsVolume(splitCandidate.mSampledFrontObjects) / viewSpaceVol;
715        splitCandidate.mVolumeBackViewCells = EvalViewCellsVolume(splitCandidate.mSampledBackObjects) / viewSpaceVol;
716
717        splitCandidate.mNumFrontViewCells = CountViewCells(splitCandidate.mSampledFrontObjects);
718        splitCandidate.mNumBackViewCells = CountViewCells(splitCandidate.mSampledBackObjects);
719
720        splitCandidate.mCorrectedFrontVolume =
721                mHierarchyManager->EvalCorrectedPvs(splitCandidate.mVolumeFrontViewCells, parentVol, avgRaysPerObject);
722       
723        splitCandidate.mCorrectedBackVolume =
724                mHierarchyManager->EvalCorrectedPvs(splitCandidate.mVolumeBackViewCells, parentVol, avgRaysPerObject);
725       
726        const float relfrontCost = splitCandidate.mCorrectedFrontVolume *
727                EvalAbsCost(splitCandidate.mFrontObjects);
728        const float relBackCost =  splitCandidate.mCorrectedBackVolume *
729                EvalAbsCost(splitCandidate.mBackObjects);
730        const float relParentCost = parentVol * EvalAbsCost(leaf->mObjects);
731
732        // compute global decrease in render cost
733        const float newRenderCost = relfrontCost + relBackCost;
734        const float renderCostDecr = relParentCost - newRenderCost;
735       
736        splitCandidate.SetRenderCostDecrease(renderCostDecr);
737
738        // increase in pvs entries
739        const int pvsEntriesIncr = EvalPvsEntriesIncr(splitCandidate,
740                                                                                                  avgRaysPerObject,
741                                                                                                  (int)pvs,
742                                                                                                  splitCandidate.mNumFrontViewCells,
743                                                                                                  splitCandidate.mNumBackViewCells);
744
745        splitCandidate.SetPvsEntriesIncr(pvsEntriesIncr);
746
747        if (0)
748        {
749                cout << "bvh volume cost"
750                 << " avg rays per object: " << avgRaysPerObject << " ratio: " << oldRatio
751                 << " parent: " << parentVol << " old vol: " << oldVolume
752                 << " frontvol: " << splitCandidate.mVolumeFrontViewCells << " corr. " << splitCandidate.mCorrectedFrontVolume
753                 << " backvol: " << splitCandidate.mVolumeBackViewCells << " corr. " << splitCandidate.mCorrectedBackVolume << endl;
754        }
755
756#ifdef GTP_DEBUG
757        Debug << "old render cost: " << oldRenderCost << endl;
758        Debug << "new render cost: " << newRenderCost << endl;
759        Debug << "render cost decrease: " << renderCostDecr << endl;
760#endif
761
762        float priority = EvalPriority(splitCandidate,
763                                                                  renderCostDecr,
764                                                                  relParentCost);
765
766        // compute global decrease in render cost
767        splitCandidate.SetPriority(priority);
768
769#if STORE_VIEWCELLS_WITH_BVH
770        if (preprocessViewCells)
771                ReleaseViewCells(*splitCandidate.mParentData.mSampledObjects);
772#endif
773
774        mEvalTimer.Exit();
775}
776
777
778int BvHierarchy::EvalPvsEntriesIncr(BvhSubdivisionCandidate &splitCandidate,
779                                                                        const float avgRaysPerObjects,
780                                                                        const int numParentViewCells,
781                                                                        const int numFrontViewCells,
782                                                                        const int numBackViewCells) //const
783{
784        const float oldPvsSize = (float)numParentViewCells;//CountViewCells(*splitCandidate.mParentData.mSampledObjects);
785        const float oldPvsRatio = (splitCandidate.mParentData.mPvs > 0) ? oldPvsSize / splitCandidate.mParentData.mPvs : 1;
786
787        const float parentPvs = splitCandidate.mParentData.mCorrectedPvs * oldPvsRatio;
788
789        const int frontViewCells = numFrontViewCells;//CountViewCells(splitCandidate.mSampledFrontObjects);
790        const int backViewCells = numBackViewCells;//CountViewCells(splitCandidate.mSampledBackObjects);
791       
792        splitCandidate.mCorrectedFrontPvs =
793                mHierarchyManager->EvalCorrectedPvs((float)frontViewCells, parentPvs, avgRaysPerObjects);
794        splitCandidate.mCorrectedBackPvs =
795                mHierarchyManager->EvalCorrectedPvs((float)backViewCells, parentPvs, avgRaysPerObjects);
796
797        if (0)
798        cout << "bvh pvs"
799                 << " avg ray contri: " << avgRaysPerObjects << " ratio: " << oldPvsRatio
800                 << " parent: " << parentPvs << " " << " old vol: " << oldPvsSize
801                 << " frontpvs: " << frontViewCells << " corr. " << splitCandidate.mCorrectedFrontPvs
802                 << " backpvs: " << frontViewCells << " corr. " << splitCandidate.mCorrectedBackPvs << endl;
803
804        return (int)(splitCandidate.mCorrectedFrontPvs + splitCandidate.mCorrectedBackPvs - parentPvs);
805}
806
807
808inline bool BvHierarchy::LocalTerminationCriteriaMet(const BvhTraversalData &tData) const
809{
810        const bool terminationCriteriaMet =
811                        (0
812                        || ((int)tData.mNode->mObjects.size() <= 1)//mTermMinObjects)
813                        //|| (data.mProbability <= mTermMinProbability)
814                        //|| (data.mNumRays <= mTermMinRays)
815                 );
816
817#ifdef _DEBUG
818        if (terminationCriteriaMet)
819        {
820                cout << "bvh local termination criteria met:" << endl;
821                cout << "objects: " << (int)tData.mNode->mObjects.size() << " (" << mTermMinObjects << ")" << endl;
822        }
823#endif
824        return terminationCriteriaMet;
825}
826
827
828inline bool BvHierarchy::GlobalTerminationCriteriaMet(const BvhTraversalData &data) const
829{
830        // note: tracking for global cost termination
831        // does not make much sense for interleaved vsp / osp partition
832        // as it is the responsibility of the hierarchy manager
833
834        const bool terminationCriteriaMet =
835                (0
836                || (mBvhStats.Leaves() >= mTermMaxLeaves)
837                //|| (mBvhStats.mGlobalCostMisses >= mTermGlobalCostMissTolerance)
838                //|| mOutOfMemory
839                );
840
841#ifdef GTP_DEBUG
842        if (terminationCriteriaMet)
843        {
844                Debug << "bvh global termination criteria met:" << endl;
845                Debug << "cost misses: " << mBvhStats.mGlobalCostMisses << " " << mTermGlobalCostMissTolerance << endl;
846                Debug << "leaves: " << mBvhStats.Leaves() << " " << mTermMaxLeaves << endl;
847        }
848#endif
849        return terminationCriteriaMet;
850}
851
852
853void BvHierarchy::EvaluateLeafStats(const BvhTraversalData &data)
854{
855        // the node became a leaf -> evaluate stats for leafs
856        BvhLeaf *leaf = data.mNode;
857       
858        ++ mCreatedLeaves;
859
860       
861        /*if (data.mProbability <= mTermMinProbability)
862        {
863                ++ mBvhStats.minProbabilityNodes;
864        }*/
865
866        ////////////////////////////////////////////
867        // depth related stuff
868
869        if (data.mDepth < mBvhStats.minDepth)
870        {
871                mBvhStats.minDepth = data.mDepth;
872        }
873
874        if (data.mDepth >= mTermMaxDepth)
875        {
876        ++ mBvhStats.maxDepthNodes;
877        }
878
879        // accumulate depth to compute average depth
880        mBvhStats.accumDepth += data.mDepth;
881
882
883        ////////////////////////////////////////////
884        // objects related stuff
885
886        // note: the sum should alwaysbe total number of objects for bvh
887        mBvhStats.objectRefs += (int)leaf->mObjects.size();
888
889        if ((int)leaf->mObjects.size() <= mTermMinObjects)
890        {
891             ++ mBvhStats.minObjectsNodes;
892        }
893
894        if (leaf->mObjects.empty())
895        {
896                ++ mBvhStats.emptyNodes;
897        }
898
899        if ((int)leaf->mObjects.size() > mBvhStats.maxObjectRefs)
900        {
901                mBvhStats.maxObjectRefs = (int)leaf->mObjects.size();
902        }
903
904        if ((int)leaf->mObjects.size() < mBvhStats.minObjectRefs)
905        {
906                mBvhStats.minObjectRefs = (int)leaf->mObjects.size();
907        }
908
909        ////////////////////////////////////////////
910        // ray related stuff
911
912        // note: this number should always accumulate to the total number of rays
913        mBvhStats.rayRefs += data.mNumRays;
914       
915        if (data.mNumRays <= mTermMinRays)
916        {
917             ++ mBvhStats.minRaysNodes;
918        }
919
920        if (data.mNumRays > mBvhStats.maxRayRefs)
921        {
922                mBvhStats.maxRayRefs = data.mNumRays;
923        }
924
925        if (data.mNumRays < mBvhStats.minRayRefs)
926        {
927                mBvhStats.minRayRefs = data.mNumRays;
928        }
929
930#ifdef _DEBUG
931        cout << "depth: " << data.mDepth << " objects: " << (int)leaf->mObjects.size()
932                 << " rays: " << data.mNumRays << " rays / objects "
933                 << (float)data.mNumRays / (float)leaf->mObjects.size() << endl;
934#endif
935}
936
937
938#if 1
939
940/// compute object boundaries using spatial mid split
941float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
942                                                                                        const int axis,
943                                                                                        ObjectContainer &objectsFront,
944                                                                                        ObjectContainer &objectsBack)
945{
946        AxisAlignedBox3 parentBox = tData.mNode->GetBoundingBox();
947
948        const float maxBox = parentBox.Max(axis);
949        const float minBox = parentBox.Min(axis);
950
951        float midPoint = (maxBox + minBox) * 0.5f;
952
953        ObjectContainer::const_iterator oit, oit_end = tData.mNode->mObjects.end();
954       
955        for (oit = tData.mNode->mObjects.begin(); oit != oit_end; ++ oit)
956        {
957                Intersectable *obj = *oit;
958                const AxisAlignedBox3 box = obj->GetBox();
959
960                const float objMid = (box.Max(axis) + box.Min(axis)) * 0.5f;
961
962                // object mailed => belongs to back objects
963                if (objMid < midPoint)
964                {
965                        objectsBack.push_back(obj);
966                }
967                else
968                {
969                        objectsFront.push_back(obj);
970                }
971        }
972
973        AxisAlignedBox3 fbox = EvalBoundingBox(objectsFront, &parentBox);
974        AxisAlignedBox3 bbox = EvalBoundingBox(objectsBack, &parentBox);
975
976        const float oldRenderCost = (float)tData.mNode->mObjects.size() * parentBox.SurfaceArea();
977        const float newRenderCost = (float)objectsFront.size() * fbox.SurfaceArea() +  (float)objectsBack.size() * bbox.SurfaceArea();
978
979        const float ratio = newRenderCost / oldRenderCost;
980        return ratio;
981}
982
983#else
984
985/// compute object partition by getting balanced objects on the left and right side
986float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
987                                                                                        const int axis,
988                                                                                        ObjectContainer &objectsFront,
989                                                                                        ObjectContainer &objectsBack)
990{
991        PrepareLocalSubdivisionCandidates(tData, axis);
992       
993        SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();
994
995        int i = 0;
996        const int border = (int)tData.mNode->mObjects.size() / 2;
997
998    for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ i)
999        {
1000                Intersectable *obj = (*cit).mObject;
1001
1002                // object mailed => belongs to back objects
1003                if (i < border)
1004                {
1005                        objectsBack.push_back(obj);
1006                }
1007                else
1008                {
1009                        objectsFront.push_back(obj);
1010                }
1011        }
1012
1013#if 1
1014        // hack: always take driving axis
1015        const float cost = (tData.mNode->GetBoundingBox().Size().DrivingAxis() == axis) ? -1.0f : 0.0f;
1016#else
1017        const float oldRenderCost = EvalAbsCost(tData.mLeaf->mObjects) / EvalProbability(tData.mSampledObjects);
1018        const float newRenderCost = EvalRenderCost(objectsFront) + EvalRenderCost(objectsBack);
1019
1020        const float cost = newRenderCost / oldRenderCost;
1021#endif
1022
1023        return cost;
1024}
1025#endif
1026
1027
1028float BvHierarchy::EvalSah(const BvhTraversalData &tData,
1029                                                   const int axis,
1030                                                   ObjectContainer &objectsFront,
1031                                                   ObjectContainer &objectsBack)
1032{
1033        // go through the lists, count the number of objects left and right
1034        // and evaluate the following cost funcion:
1035        // C = ct_div_ci  + (ol + or) / queries
1036        PrepareLocalSubdivisionCandidates(tData, axis);
1037
1038        const float totalRenderCost = EvalAbsCost(tData.mNode->mObjects);
1039        float objectsLeft = 0, objectsRight = totalRenderCost;
1040 
1041        const AxisAlignedBox3 nodeBbox = tData.mNode->GetBoundingBox();
1042        const float boxArea = nodeBbox.SurfaceArea();
1043
1044        float minSum = 1e20f;
1045 
1046        float minBorder = nodeBbox.Max(axis);
1047        float maxBorder = nodeBbox.Min(axis);
1048
1049        float areaLeft = 0, areaRight = 0;
1050
1051        SortableEntryContainer::const_iterator currentPos =
1052                mSubdivisionCandidates->begin();
1053       
1054        vector<float> bordersRight;
1055
1056        // we keep track of both borders of the bounding boxes =>
1057        // store the events in descending order
1058
1059        bordersRight.resize(mSubdivisionCandidates->size());
1060
1061        SortableEntryContainer::reverse_iterator rcit =
1062                mSubdivisionCandidates->rbegin(), rcit_end = mSubdivisionCandidates->rend();
1063
1064        vector<float>::reverse_iterator rbit = bordersRight.rbegin();
1065
1066        for (; rcit != rcit_end; ++ rcit, ++ rbit)
1067        {
1068                Intersectable *obj = (*rcit).mObject;
1069                const AxisAlignedBox3 obox = obj->GetBox();
1070
1071                if (obox.Min(axis) < minBorder)
1072                {
1073                        minBorder = obox.Min(axis);
1074                }
1075
1076                (*rbit) = minBorder;
1077        }
1078
1079        // record surface areas during the sweep
1080        float al = 0;
1081        float ar = boxArea;
1082
1083        vector<float>::const_iterator bit = bordersRight.begin();
1084        SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();
1085
1086        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit)
1087        {
1088                Intersectable *obj = (*cit).mObject;
1089
1090#if USE_BETTER_RENDERCOST_EST
1091                const float renderCost = ViewCellsManager::EvalRenderCost(obj);
1092               
1093                objectsLeft += renderCost;
1094                objectsRight -= renderCost;
1095
1096                const bool noValidSplit = ((objectsLeft <= Limits::Small) || (objectsRight <= Limits::Small));
1097
1098#else
1099                ++ objectsLeft;
1100                -- objectsRight;
1101
1102                const bool noValidSplit = !objectsLeft || !objectsRight;
1103#endif
1104                const AxisAlignedBox3 obox = obj->GetBox();
1105
1106                // the borders of the bounding boxes have changed
1107                if (obox.Max(axis) > maxBorder)
1108                {
1109                        maxBorder = obox.Max(axis);
1110                }
1111
1112                minBorder = (*bit);
1113
1114                AxisAlignedBox3 lbox = nodeBbox;
1115                AxisAlignedBox3 rbox = nodeBbox;
1116
1117                lbox.SetMax(axis, maxBorder);
1118                rbox.SetMin(axis, minBorder);
1119
1120                al = lbox.SurfaceArea();
1121                ar = rbox.SurfaceArea();
1122
1123                // should use classical approach here ...
1124#if BOUND_RENDERCOST
1125                const float rcLeft = std::max(objectsLeft, MIN_RENDERCOST);
1126                const float rcRight = std::max(objectsRight, MIN_RENDERCOST);
1127
1128                const float sum = noValidSplit ? 1e25f : objectsLeft * al + objectsRight * ar;
1129#else
1130
1131                const float sum = noValidSplit ? 1e25f : objectsLeft * al + objectsRight * ar;
1132#endif
1133                /*cout << "pos=" << (*cit).mPos << "\t q=(" << objectsLeft << "," << objectsRight <<")\t r=("
1134                         << lbox.SurfaceArea() << "," << rbox.SurfaceArea() << ")" << endl;
1135                cout << "minborder: " << minBorder << " maxborder: " << maxBorder << endl;
1136            cout << "cost= " << sum << endl;*/
1137       
1138                if (sum < minSum)
1139                {       
1140                        minSum = sum;
1141                        areaLeft = al;
1142                        areaRight = ar;
1143
1144                        // objects belong to left side now
1145                        for (; currentPos != (cit + 1); ++ currentPos);
1146                }
1147        }
1148
1149        ////////////
1150        //-- assign object to front and back volume
1151
1152        // belongs to back bv
1153        for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
1154                objectsBack.push_back((*cit).mObject);
1155       
1156        // belongs to front bv
1157        for (cit = currentPos; cit != cit_end; ++ cit)
1158                objectsFront.push_back((*cit).mObject);
1159
1160        float newCost = minSum / boxArea;
1161        float ratio = newCost / totalRenderCost;
1162 
1163#ifdef GTP_DEBUG
1164        cout << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of "
1165                 << (int)tData.mNode->mObjects.size() << ")\t area=("
1166                 << areaLeft << ", " << areaRight << ", " << boxArea << ")" << endl
1167                 << "cost= " << newCost << " oldCost=" << totalRenderCost / boxArea << endl;
1168#endif
1169
1170        return ratio;
1171}
1172
1173
1174
1175float BvHierarchy::EvalSahWithTigherBbox(const BvhTraversalData &tData,
1176                                                                                 const int axis,
1177                                                                                 ObjectContainer &objectsFront,
1178                                                                                 ObjectContainer &objectsBack)
1179{
1180        // go through the lists, count the number of objects left and right
1181        // and evaluate the following cost funcion:
1182        // C = ct_div_ci  + (ol + or) / queries
1183        PrepareLocalSubdivisionCandidates(tData, axis);
1184
1185        const float totalRenderCost = EvalAbsCost(tData.mNode->mObjects);
1186        float objectsLeft = 0, objectsRight = totalRenderCost;
1187 
1188        const AxisAlignedBox3 nodeBbox = tData.mNode->GetBoundingBox();
1189
1190        const float minBox = nodeBbox.Min(axis);
1191        const float maxBox = nodeBbox.Max(axis);
1192        const float boxArea = nodeBbox.SurfaceArea();
1193
1194        float minSum = 1e20f;
1195 
1196        Vector3 minBorder = nodeBbox.Max();
1197        Vector3 maxBorder = nodeBbox.Min();
1198
1199        float areaLeft = 0, areaRight = 0;
1200
1201        SortableEntryContainer::const_iterator currentPos =
1202                mSubdivisionCandidates->begin();
1203       
1204        vector<Vector3> bordersRight;
1205
1206        // we keep track of both borders of the bounding boxes =>
1207        // store the events in descending order
1208        bordersRight.resize(mSubdivisionCandidates->size());
1209
1210        SortableEntryContainer::reverse_iterator rcit =
1211                mSubdivisionCandidates->rbegin(), rcit_end =
1212                mSubdivisionCandidates->rend();
1213
1214        vector<Vector3>::reverse_iterator rbit = bordersRight.rbegin();
1215
1216        for (; rcit != rcit_end; ++ rcit, ++ rbit)
1217        {
1218                Intersectable *obj = (*rcit).mObject;
1219                const AxisAlignedBox3 obox = obj->GetBox();
1220
1221                for (int i = 0; i < 3; ++ i)
1222                {
1223                        if (obox.Min(i) < minBorder[i])
1224                        {
1225                                minBorder[i] = obox.Min(i);
1226                        }
1227                }
1228
1229                (*rbit) = minBorder;
1230        }
1231
1232        // temporary surface areas
1233        float al = 0;
1234        float ar = boxArea;
1235
1236        vector<Vector3>::const_iterator bit = bordersRight.begin();
1237        SortableEntryContainer::const_iterator cit, cit_end =
1238                mSubdivisionCandidates->end();
1239
1240        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit)
1241        {
1242                Intersectable *obj = (*cit).mObject;
1243
1244                const float renderCost = ViewCellsManager::EvalRenderCost(obj);
1245               
1246                objectsLeft += renderCost;
1247                objectsRight -= renderCost;
1248
1249                const AxisAlignedBox3 obox = obj->GetBox();
1250
1251                AxisAlignedBox3 lbox = nodeBbox;
1252                AxisAlignedBox3 rbox = nodeBbox;
1253       
1254                // the borders of the left bounding box have changed
1255                for (int i = 0; i < 3; ++ i)
1256                {
1257                        if (obox.Max(i) > maxBorder[i])
1258                        {
1259                                maxBorder[i] = obox.Max(i);
1260                        }
1261                }
1262
1263                minBorder = (*bit);
1264
1265                lbox.SetMax(maxBorder);
1266                rbox.SetMin(minBorder);
1267
1268                al = lbox.SurfaceArea();
1269                ar = rbox.SurfaceArea();
1270       
1271                const bool noValidSplit = ((objectsLeft <= Limits::Small) || (objectsRight <= Limits::Small));
1272                const float sum =  noValidSplit ? 1e25 : objectsLeft * al + objectsRight * ar;
1273     
1274                /*cout << "pos=" << (*cit).mPos << "\t q=(" << objectsLeft << "," << objectsRight <<")\t r=("
1275                         << lbox.SurfaceArea() << "," << rbox.SurfaceArea() << ")" << endl;
1276                cout << "minborder: " << minBorder << " maxborder: " << maxBorder << endl;
1277            cout << "cost= " << sum << endl;*/
1278       
1279                if (sum < minSum)
1280                {       
1281                        minSum = sum;
1282                        areaLeft = al;
1283                        areaRight = ar;
1284
1285                        // objects belong to left side now
1286                        for (; currentPos != (cit + 1); ++ currentPos);
1287                }
1288        }
1289
1290        /////////////
1291        //-- assign object to front and back volume
1292
1293        // belongs to back bv
1294        for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
1295                objectsBack.push_back((*cit).mObject);
1296       
1297        // belongs to front bv
1298        for (cit = currentPos; cit != cit_end; ++ cit)
1299                objectsFront.push_back((*cit).mObject);
1300
1301        float newCost = minSum / boxArea;
1302        float ratio = newCost / totalRenderCost;
1303 
1304#ifdef GTP_DEBUG
1305        cout << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of "
1306                 << (int)tData.mNode->mObjects.size() << ")\t area=("
1307                 << areaLeft << ", " << areaRight << ", " << boxArea << ")" << endl
1308                 << "cost= " << newCost << " oldCost=" << totalRenderCost / boxArea << endl;
1309#endif
1310
1311        return ratio;
1312}
1313
1314
1315static bool PrepareOutput(const int axis,
1316                                                  const int leaves,
1317                                                  ofstream &sumStats,
1318                                                  ofstream &vollStats,
1319                                                  ofstream &volrStats)
1320{
1321        if ((axis == 0) && (leaves > 0) && (leaves < 90))
1322        {
1323                char str[64];   
1324                sprintf(str, "tmp/bvh_heur_sum-%04d.log", leaves);
1325                sumStats.open(str);
1326                sprintf(str, "tmp/bvh_heur_voll-%04d.log", leaves);
1327                vollStats.open(str);
1328                sprintf(str, "tmp/bvh_heur_volr-%04d.log", leaves);
1329                volrStats.open(str);
1330        }
1331
1332        return sumStats.is_open() && vollStats.is_open() && volrStats.is_open();
1333}
1334
1335
1336static void PrintHeuristics(const float objectsRight,
1337                                                        const float sum,
1338                                                        const float volLeft,
1339                                                        const float volRight,
1340                                                        const float viewSpaceVol,
1341                                                        ofstream &sumStats,
1342                                                        ofstream &vollStats,
1343                                                        ofstream &volrStats)
1344{
1345        sumStats
1346                << "#Position\n" << objectsRight << endl
1347                << "#Sum\n" << sum / viewSpaceVol << endl
1348                << "#Vol\n" << (volLeft +  volRight) / viewSpaceVol << endl;
1349
1350        vollStats
1351                << "#Position\n" << objectsRight << endl
1352                << "#Vol\n" << volLeft / viewSpaceVol << endl;
1353
1354        volrStats
1355                << "#Position\n" << objectsRight << endl
1356                << "#Vol\n" << volRight / viewSpaceVol << endl;
1357}
1358
1359
1360float BvHierarchy::EvalLocalCostHeuristics(const BvhTraversalData &tData,
1361                                                                                   const int axis,
1362                                                                                   ObjectContainer &objectsFront,
1363                                                                                   ObjectContainer &objectsBack)
1364{
1365        /////////////////////////////////////////////
1366        //-- go through the lists, count the number of objects
1367        //-- left and right and evaluate the cost funcion
1368
1369        // prepare the heuristics by setting mailboxes and counters
1370        const float totalVol = PrepareHeuristics(tData, axis);
1371       
1372        // local helper variables
1373        float volLeft = 0;
1374        float volRight = totalVol;
1375       
1376        const float nTotalObjects = EvalAbsCost(tData.mNode->mObjects);
1377        float nObjectsLeft = 0;
1378        float nObjectsRight = nTotalObjects;
1379
1380        const float viewSpaceVol =
1381                mViewCellsManager->GetViewSpaceBox().GetVolume();
1382
1383        SortableEntryContainer::const_iterator backObjectsStart =
1384                mSubdivisionCandidates->begin();
1385
1386        /////////////////////////////////
1387        //-- the parameters for the current optimum
1388
1389        float volBack = volLeft;
1390        float volFront = volRight;
1391        float newRenderCost = nTotalObjects * totalVol;
1392
1393#ifdef GTP_DEBUG
1394        ofstream sumStats;
1395        ofstream vollStats;
1396        ofstream volrStats;
1397
1398        const bool printStats = PrepareOutput(axis,
1399                                                                                  mBvhStats.Leaves(),
1400                                                                                  sumStats,
1401                                                                                  vollStats,
1402                                                                                  volrStats);
1403#endif
1404
1405        ///////////////////////
1406        //-- the sweep heuristics
1407        //-- traverse through events and find best split plane
1408
1409        SortableEntryContainer::const_iterator cit,
1410                cit_end = cit_end = mSubdivisionCandidates->end();
1411
1412        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit)
1413        {
1414                Intersectable *object = (*cit).mObject;
1415       
1416                // evaluate change in l and r volume
1417                // voll = view cells that see only left node (i.e., left pvs)
1418                // volr = view cells that see only right node (i.e., right pvs)
1419                EvalHeuristicsContribution(object, volLeft, volRight);
1420
1421#if USE_BETTER_RENDERCOST_EST
1422
1423                const float rc = ViewCellsManager::EvalRenderCost(object);
1424               
1425                nObjectsLeft += rc;
1426                nObjectsRight -= rc;
1427       
1428#else
1429                ++ nObjectsLeft;
1430                -- nObjectsRight;
1431#endif
1432       
1433                // split is only valid if #objects on left and right is not zero
1434                const bool noValidSplit = (nObjectsRight <= Limits::Small);
1435
1436#if BOUND_RENDERCOST
1437
1438                const float rcLeft = max(nObjectsLeft, MIN_RENDERCOST);
1439                const float rcRight = max(nObjectsRight, MIN_RENDERCOST);
1440               
1441                // the heuristics
1442            const float sum = noValidSplit ?
1443                        1e25f : volLeft * (float)rcLeft + volRight * (float)rcRight;
1444#else
1445       
1446                // the heuristics
1447            const float sum = noValidSplit ?
1448                        1e25f : volLeft * (float)nObjectsLeft + volRight * (float)nObjectsRight;
1449#endif
1450               
1451#ifdef GTP_DEBUG
1452                if (printStats)
1453                {
1454                        PrintHeuristics(nObjectsRight, sum, volLeft,
1455                                                        volRight, viewSpaceVol,
1456                                                        sumStats, vollStats, volrStats);
1457                }
1458#endif
1459
1460                if (sum < newRenderCost)
1461                {
1462                        newRenderCost = sum;
1463
1464                        volBack = volLeft;
1465                        volFront = volRight;
1466
1467                        // objects belongs to left side now
1468                        for (; backObjectsStart != (cit + 1); ++ backObjectsStart);
1469                }
1470        }
1471
1472        ////////////////////////////////////////
1473        //-- assign object to front and back volume
1474
1475        // belongs to back bv
1476        for (cit = mSubdivisionCandidates->begin(); cit != backObjectsStart; ++ cit)
1477        {
1478                objectsBack.push_back((*cit).mObject);
1479        }
1480        // belongs to front bv
1481        for (cit = backObjectsStart; cit != cit_end; ++ cit)
1482        {
1483                objectsFront.push_back((*cit).mObject);
1484        }
1485
1486        // render cost of the old parent
1487        const float oldRenderCost = (float)nTotalObjects * totalVol + Limits::Small;
1488        // the relative cost ratio
1489        const float ratio = newRenderCost / oldRenderCost;
1490
1491#ifdef GTP_DEBUG
1492        Debug << "\n§§§§ bvh eval const decrease §§§§" << endl
1493                  << "back pvs: " << (int)objectsBack.size() << " front pvs: "
1494                  << (int)objectsFront.size() << " total pvs: " << nTotalObjects << endl
1495                  << "back p: " << volBack / viewSpaceVol << " front p "
1496                  << volFront / viewSpaceVol << " p: " << totalVol / viewSpaceVol << endl
1497                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: "
1498                  << newRenderCost / viewSpaceVol << endl
1499                  << "render cost decrease: "
1500                  << oldRenderCost / viewSpaceVol - newRenderCost / viewSpaceVol << endl;
1501#endif
1502
1503        return ratio;
1504}
1505
1506
1507void BvHierarchy::PrepareLocalSubdivisionCandidates(const BvhTraversalData &tData,
1508                                                                                                        const int axis)                                                                                 
1509{
1510        mSortTimer.Entry();
1511       
1512        //-- insert object queries
1513        ObjectContainer *objects = mUseGlobalSorting ?
1514                tData.mSortedObjects[axis] : &tData.mNode->mObjects;
1515
1516        CreateLocalSubdivisionCandidates(*objects, &mSubdivisionCandidates, !mUseGlobalSorting, axis);
1517       
1518        mSortTimer.Exit();
1519}
1520
1521
1522void BvHierarchy::CreateLocalSubdivisionCandidates(const ObjectContainer &objects,
1523                                                                                                  SortableEntryContainer **subdivisionCandidates,
1524                                                                                                  const bool sortEntries,
1525                                                                                                  const int axis)
1526{
1527        (*subdivisionCandidates)->clear();
1528
1529        // compute requested size and look if subdivision candidate has to be recomputed
1530        const int requestedSize = (int)objects.size();
1531       
1532        // creates a sorted split candidates array
1533        if ((*subdivisionCandidates)->capacity() > 500000 &&
1534                requestedSize < (int)((*subdivisionCandidates)->capacity() / 10) )
1535        {
1536        delete (*subdivisionCandidates);
1537                (*subdivisionCandidates) = new SortableEntryContainer;
1538        }
1539
1540        (*subdivisionCandidates)->reserve(requestedSize);
1541
1542        ObjectContainer::const_iterator oit, oit_end = objects.end();
1543
1544        for (oit = objects.begin(); oit < oit_end; ++ oit)
1545        {
1546                (*subdivisionCandidates)->push_back(SortableEntry(*oit, (*oit)->GetBox().Center(axis)));
1547        }
1548
1549        if (sortEntries)
1550        {       // no presorted candidate list
1551                stable_sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
1552                //sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
1553        }
1554}
1555
1556
1557const BvhStatistics &BvHierarchy::GetStatistics() const
1558{
1559        return mBvhStats;
1560}
1561
1562
1563float BvHierarchy::PrepareHeuristics(const BvhTraversalData &tData,
1564                                                                         const int axis)
1565{       
1566        BvhLeaf *leaf = tData.mNode;
1567        float vol = 0;
1568
1569    // sort so we can use a sweep from right to left
1570        PrepareLocalSubdivisionCandidates(tData, axis);
1571       
1572        // collect and mark the view cells as belonging to front pvs
1573        ViewCellContainer viewCells;
1574
1575        const bool setCounter = true;
1576        const bool onlyUnmailed = true;
1577
1578       
1579        CollectViewCells(*tData.mSampledObjects,
1580                                         viewCells,
1581                                         setCounter,
1582                                         onlyUnmailed);
1583
1584        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
1585
1586        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
1587        {
1588#if USE_VOLUMES_FOR_HEURISTICS
1589                const float volIncr = (*vit)->GetVolume();
1590#else
1591                const float volIncr = 1.0f;
1592#endif
1593                vol += volIncr;
1594        }
1595
1596        // mail view cells that go from front node to back node
1597        ViewCell::NewMail();
1598       
1599        return vol;
1600}
1601
1602///////////////////////////////////////////////////////////
1603
1604
1605void BvHierarchy::EvalHeuristicsContribution(Intersectable *obj,
1606                                                                                         float &volLeft,
1607                                                                                         float &volRight)
1608{
1609        // collect all view cells associated with this objects
1610        // (also multiple times, if they are pierced by several rays)
1611        ViewCellContainer viewCells;
1612
1613        const bool useMailboxing = false;
1614        const bool setCounter = false;
1615        const bool onlyUnmailedRays = true;
1616
1617        CollectViewCells(obj, viewCells, useMailboxing, setCounter, onlyUnmailedRays);
1618
1619        // classify view cells and compute volume contri accordingly
1620        // possible view cell classifications:
1621        // view cell mailed => view cell can be seen from left child node
1622        // view cell counter > 0 view cell can be seen from right child node
1623        // combined: view cell volume belongs to both nodes
1624        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
1625       
1626        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
1627        {
1628                // view cells can also be seen from left child node
1629                ViewCell *viewCell = *vit;
1630
1631#if USE_VOLUMES_FOR_HEURISTICS
1632                const float vol = viewCell->GetVolume();
1633#else
1634                const float vol = 1.0f;
1635#endif
1636                if (!viewCell->Mailed())
1637                {
1638                        viewCell->Mail();
1639                        // we now see view cell from both nodes
1640                        // => add volume to left node
1641                        volLeft += vol;
1642                }
1643
1644                // last reference into the right node
1645                if (-- viewCell->mCounter == 0)
1646                {       
1647                        // view cell was previously seen from both nodes  =>
1648                        // remove volume from right node
1649                        volRight -= vol;
1650                }
1651        }
1652}
1653
1654
1655void BvHierarchy::SetViewCellsManager(ViewCellsManager *vcm)
1656{
1657        mViewCellsManager = vcm;
1658}
1659
1660
1661AxisAlignedBox3 BvHierarchy::GetBoundingBox() const
1662{
1663        return mBoundingBox;
1664}
1665
1666
1667float BvHierarchy::SelectObjectPartition(const BvhTraversalData &tData,
1668                                                                                 ObjectContainer &frontObjects,
1669                                                                                 ObjectContainer &backObjects,
1670                                                                                 bool useVisibilityBasedHeuristics)
1671{
1672        mSplitTimer.Entry();
1673
1674        if (mIsInitialSubdivision)
1675        {
1676                ApplyInitialSplit(tData, frontObjects, backObjects);
1677                return 0;
1678        }
1679
1680        ObjectContainer nFrontObjects[3];
1681        ObjectContainer nBackObjects[3];
1682        float nCostRatio[3];
1683
1684        int sAxis = 0;
1685        int bestAxis = -1;
1686
1687        if (mOnlyDrivingAxis)
1688        {
1689                const AxisAlignedBox3 box = tData.mNode->GetBoundingBox();
1690                sAxis = box.Size().DrivingAxis();
1691        }
1692
1693        // if #rays high consider only use a subset of the rays for
1694        // visibility based heuristics
1695        VssRay::NewMail();
1696
1697        /*if ((mMaxTests < tData.mNumRays) && mUseCostHeuristics && useVisibilityBasedHeuristics)
1698        {
1699                VssRayContainer rays;
1700
1701                // maximal 2 objects share the same ray
1702                rays.reserve(tData.mNumRays * 2);
1703                CollectRays(tData.mNode->mObjects, rays);
1704
1705                const float prop = (float)mMaxTests / (float)rays.size();
1706
1707                VssRayContainer::const_iterator rit, rit_end = rays.end();
1708
1709                // mail rays which will not be considered
1710                for (rit = rays.begin(); rit != rit_end; ++ rit)
1711                {
1712                        if (Random(1.0f) > prop)
1713                        {
1714                                (*rit)->Mail();
1715                        }
1716                }               
1717        }*/
1718
1719        ////////////////////////////////////
1720        //-- evaluate split cost for all three axis
1721       
1722        for (int axis = 0; axis < 3; ++ axis)
1723        {
1724                if (!mOnlyDrivingAxis || (axis == sAxis))
1725                {
1726                        if (mUseCostHeuristics)
1727                        {
1728                                //////////////////////////////////
1729                //-- split objects using heuristics
1730                               
1731                                if (useVisibilityBasedHeuristics)
1732                                {
1733                                        ///////////
1734                                        //-- heuristics using objects weighted by view cells volume
1735                                        nCostRatio[axis] =
1736                                                EvalLocalCostHeuristics(tData,
1737                                                                                                axis,
1738                                                                                                nFrontObjects[axis],
1739                                                                                                nBackObjects[axis]);
1740                                }
1741                                else
1742                                {       
1743                                        //////////////////
1744                                        //-- view cells not constructed yet     => use surface area heuristic                   
1745                                        nCostRatio[axis] = EvalSah(tData,
1746                                                                                           axis,
1747                                                                                           nFrontObjects[axis],
1748                                                                                           nBackObjects[axis]);
1749                                }
1750                        }
1751                        else
1752                        {
1753                                //-- split objects using some simple criteria
1754                                nCostRatio[axis] =
1755                                        EvalLocalObjectPartition(tData, axis, nFrontObjects[axis], nBackObjects[axis]);
1756                        }
1757
1758                        // avoid splits in degenerate axis with high penalty
1759                        if (1 &&
1760                                (tData.mNode->GetBoundingBox().Size(axis) < 0.0001))//Limits::Small))
1761                        {
1762                                nCostRatio[axis] += 9999;
1763                        }
1764
1765                        if ((bestAxis == -1) || (nCostRatio[axis] < nCostRatio[bestAxis]))
1766                        {
1767                                bestAxis = axis;
1768                        }
1769                }
1770        }
1771
1772    ////////////////
1773        //-- assign values
1774
1775        frontObjects = nFrontObjects[bestAxis];
1776        backObjects = nBackObjects[bestAxis];
1777
1778        mSplitTimer.Exit();
1779
1780        //cout << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
1781        return nCostRatio[bestAxis];
1782}
1783
1784
1785int BvHierarchy::AssociateObjectsWithRays(const VssRayContainer &rays) const
1786{
1787        int nRays = 0;
1788        VssRayContainer::const_iterator rit, rit_end = rays.end();
1789
1790        VssRay *lastVssRay = NULL;
1791
1792        VssRay::NewMail();
1793
1794    for (rit = rays.begin(); rit != rays.end(); ++ rit)
1795        {
1796                VssRay *ray = (*rit);
1797
1798                // filter out double rays (last ray the same as this ray)
1799                if (
1800                        !lastVssRay ||
1801                        !(ray->mOrigin == lastVssRay->mTermination) ||
1802                        !(ray->mTermination == lastVssRay->mOrigin))
1803                {
1804                        lastVssRay = ray;
1805                        //cout << "j";
1806                        if (ray->mTerminationObject)
1807                        {
1808                                ray->mTerminationObject->GetOrCreateRays()->push_back(ray);
1809                                if (!ray->Mailed())
1810                                {
1811                                        ray->Mail();
1812                                        ++ nRays;
1813                                }
1814                        }
1815
1816#if COUNT_ORIGIN_OBJECTS
1817
1818                        if (ray->mOriginObject)
1819                        {
1820                                //cout << "o";
1821                                ray->mOriginObject->GetOrCreateRays()->push_back(ray);
1822
1823                                if (!ray->Mailed())
1824                                {
1825                                        ray->Mail();
1826                                        ++ nRays;
1827                                }
1828                        }
1829#endif
1830                }
1831
1832        }
1833
1834        return nRays;
1835}
1836
1837
1838void BvHierarchy::PrintSubdivisionStats(const SubdivisionCandidate &sc)
1839{
1840        const float costDecr = sc.GetRenderCostDecrease();     
1841
1842        mSubdivisionStats
1843                        << "#Leaves\n" << mBvhStats.Leaves() << endl
1844                        << "#RenderCostDecrease\n" << costDecr << endl
1845                        << "#TotalRenderCost\n" << mTotalCost << endl
1846                        << "#EntriesInPvs\n" << mPvsEntries << endl;
1847}
1848
1849
1850void BvHierarchy::CollectRays(const ObjectContainer &objects,
1851                                                          VssRayContainer &rays) const
1852{
1853        VssRay::NewMail();
1854        ObjectContainer::const_iterator oit, oit_end = objects.end();
1855
1856        // evaluate reverse pvs and view cell volume on left and right cell
1857        // note: should I take all leaf objects or rather the objects hit by rays?
1858        for (oit = objects.begin(); oit != oit_end; ++ oit)
1859        {
1860                Intersectable *obj = *oit;
1861                VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();
1862
1863                for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
1864                {
1865                        VssRay *ray = (*rit);
1866
1867                        if (!ray->Mailed())
1868                        {
1869                                ray->Mail();
1870                                rays.push_back(ray);
1871                        }
1872                }
1873        }
1874}
1875
1876
1877float BvHierarchy::EvalSahCost(BvhLeaf *leaf) const
1878{
1879        ////////////////
1880        //-- surface area heuristics
1881
1882        // eraly exit
1883        //if (leaf->mObjects.empty())   return 0.0f;
1884
1885        const AxisAlignedBox3 box = GetBoundingBox(leaf);
1886        const float area = box.SurfaceArea();
1887        const float viewSpaceArea = mViewCellsManager->GetViewSpaceBox().SurfaceArea();
1888
1889        return EvalAbsCost(leaf->mObjects) * area / viewSpaceArea;
1890}
1891
1892
1893float BvHierarchy::EvalRenderCost(const ObjectContainer &objects)// const
1894{       
1895        ///////////////
1896        //-- render cost heuristics
1897
1898        const float objRenderCost = EvalAbsCost(objects);
1899
1900        const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
1901
1902        // probability that view point lies in a view cell which sees this node
1903        const float p = EvalViewCellsVolume(objects) / viewSpaceVol;
1904       
1905        return objRenderCost * p;
1906}
1907
1908
1909float BvHierarchy::EvalProbability(const ObjectContainer &objects)// const
1910{       
1911        const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
1912       
1913        // probability that view point lies in a view cell which sees this node
1914        return EvalViewCellsVolume(objects) / viewSpaceVol;
1915}
1916
1917
1918AxisAlignedBox3 BvHierarchy::EvalBoundingBox(const ObjectContainer &objects,
1919                                                                                         const AxisAlignedBox3 *parentBox) const
1920{
1921        // if there are no objects in this box, box size is set to parent box size.
1922        // Question: Invalidate box instead?
1923        if (parentBox && objects.empty())
1924                return *parentBox;
1925
1926        AxisAlignedBox3 box;
1927        box.Initialize();
1928
1929        ObjectContainer::const_iterator oit, oit_end = objects.end();
1930
1931        for (oit = objects.begin(); oit != oit_end; ++ oit)
1932        {
1933                Intersectable *obj = *oit;
1934                // grow bounding box to include all objects
1935                box.Include(obj->GetBox());
1936        }
1937
1938        return box;
1939}
1940
1941
1942void BvHierarchy::CollectLeaves(BvhNode *root, vector<BvhLeaf *> &leaves) const
1943{
1944        stack<BvhNode *> nodeStack;
1945        nodeStack.push(root);
1946
1947        while (!nodeStack.empty())
1948        {
1949                BvhNode *node = nodeStack.top();
1950                nodeStack.pop();
1951
1952                if (node->IsLeaf())
1953                {
1954                        BvhLeaf *leaf = (BvhLeaf *)node;
1955                        leaves.push_back(leaf);
1956                }
1957                else
1958                {
1959                        BvhInterior *interior = (BvhInterior *)node;
1960
1961                        nodeStack.push(interior->GetBack());
1962                        nodeStack.push(interior->GetFront());
1963                }
1964        }
1965}
1966
1967
1968void BvHierarchy::CollectNodes(BvhNode *root, vector<BvhNode *> &nodes) const
1969{
1970        stack<BvhNode *> nodeStack;
1971        nodeStack.push(root);
1972
1973        while (!nodeStack.empty())
1974        {
1975                BvhNode *node = nodeStack.top();
1976                nodeStack.pop();
1977
1978                nodes.push_back(node);
1979               
1980                if (!node->IsLeaf())
1981                {
1982                        BvhInterior *interior = (BvhInterior *)node;
1983
1984                        nodeStack.push(interior->GetBack());
1985                        nodeStack.push(interior->GetFront());
1986                }
1987        }
1988}
1989
1990
1991AxisAlignedBox3 BvHierarchy::GetBoundingBox(BvhNode *node) const
1992{
1993        return node->GetBoundingBox();
1994}
1995
1996
1997int BvHierarchy::CollectViewCells(const ObjectContainer &objects,
1998                                                                  ViewCellContainer &viewCells,
1999                                                                  const bool setCounter,
2000                                                                  const bool onlyUnmailedRays)// const
2001{
2002        ViewCell::NewMail();
2003
2004        ObjectContainer::const_iterator oit, oit_end = objects.end();
2005
2006        // use mailing to avoid dublicates
2007        const bool useMailBoxing = true;
2008
2009        int numRays = 0;
2010        // loop through all object and collect view cell pvs of this node
2011        for (oit = objects.begin(); oit != oit_end; ++ oit)
2012        {
2013                // use mailing to avoid duplicates
2014                numRays += CollectViewCells(*oit, viewCells, useMailBoxing, setCounter, onlyUnmailedRays);
2015        }
2016
2017        return numRays;
2018}
2019
2020
2021#if STORE_VIEWCELLS_WITH_BVH
2022
2023
2024void BvHierarchy::ReleaseViewCells(const ObjectContainer &objects)
2025{
2026        ObjectContainer::const_iterator oit, oit_end = objects.end();
2027
2028        for (oit = objects.begin(); oit != oit_end; ++ oit)
2029        {
2030                (*oit)->DelViewCells();
2031        }
2032}
2033
2034
2035void BvHierarchy::AssociateViewCellsWithObjects(const ObjectContainer &objects) const
2036{
2037        ObjectContainer::const_iterator oit, oit_end = objects.end();
2038
2039        const bool useMailBoxing = true;
2040        VssRay::NewMail();
2041       
2042        for (oit = objects.begin(); oit != oit_end; ++ oit)
2043        {
2044                        ViewCell::NewMail();
2045                        // use mailing to avoid duplicates
2046                        AssociateViewCellsWithObject(*oit, useMailBoxing);
2047        }
2048}
2049
2050
2051int BvHierarchy::AssociateViewCellsWithObject(Intersectable *obj, const bool useMailBoxing) const
2052{
2053        int nRays = 0;
2054
2055        if (!obj->GetOrCreateViewCells()->empty())
2056        {
2057                cerr << "AssociateViewCellsWithObject: view cells cache not working" << endl;
2058        }
2059
2060        ViewCellContainer *objViewCells = obj->GetOrCreateViewCells();
2061        VssRayContainer *vssRays = obj->GetOrCreateRays();
2062
2063        VssRayContainer::const_iterator rit, rit_end = vssRays->end();
2064
2065        // fill cache
2066        for (rit = vssRays->begin(); rit < rit_end; ++ rit)
2067        {
2068                VssRay *ray = (*rit);
2069
2070                //      if (onlyUnmailedRays && ray->Mailed())
2071                //              continue;
2072                mHierarchyManager->mVspTree->GetViewCells(*ray, *objViewCells);
2073               
2074                if (!useMailBoxing || !ray->Mailed())
2075                {
2076                        if (useMailBoxing)
2077                                ray->Mail();
2078
2079                        ++ nRays;
2080                }
2081        }
2082
2083        return nRays;
2084}
2085
2086
2087
2088int BvHierarchy::CountViewCells(Intersectable *obj) //const
2089{
2090        ViewCellContainer *viewCells = obj->GetOrCreateViewCells();
2091
2092        if (obj->GetOrCreateViewCells()->empty())
2093        {
2094                //cerr << "h";//CountViewCells: view cells empty, view cells cache not working" << endl;
2095                return CountViewCellsFromRays(obj);
2096        }
2097       
2098        int result = 0;
2099
2100        ViewCellContainer::const_iterator vit, vit_end = viewCells->end();
2101       
2102        for (vit = viewCells->begin(); vit != vit_end; ++ vit)
2103        {
2104                ViewCell *vc = *vit;
2105
2106                // store view cells
2107                if (!vc->Mailed())
2108                {
2109                        vc->Mail();
2110                        ++ result;
2111                }
2112        }
2113
2114        return result;
2115}
2116
2117
2118int BvHierarchy::CollectViewCells(Intersectable *obj,
2119                                                                  ViewCellContainer &viewCells,
2120                                                                  const bool useMailBoxing,
2121                                                                  const bool setCounter,
2122                                                                  const bool onlyUnmailedRays)// const
2123{
2124        if (obj->GetOrCreateViewCells()->empty())
2125        {
2126                //cerr << "g";//CollectViewCells: view cells empty, view cells cache not working" << endl;
2127                return CollectViewCellsFromRays(obj, viewCells, useMailBoxing, setCounter, onlyUnmailedRays);
2128        }
2129
2130        mCollectTimer.Entry();
2131
2132        ViewCellContainer *objViewCells = obj->GetOrCreateViewCells();
2133
2134        ///////////
2135        //-- use view cells cache
2136
2137        // fastest way: just copy to the end
2138        //if (!useMailBoxing)   viewCells.insert(viewCells.end(), objViewCells->begin(), objViewCells->end());
2139
2140        // loop through view cells
2141        // matt: probably slow to insert  view cells one by one
2142        ViewCellContainer::const_iterator vit, vit_end = objViewCells->end();
2143
2144        for (vit = objViewCells->begin(); vit != vit_end; ++ vit)
2145        {
2146                ViewCell *vc = *vit;
2147
2148                // store view cells
2149                if (!useMailBoxing || !vc->Mailed())
2150                {
2151                        if (useMailBoxing)
2152                        {
2153                                // view cell not mailed
2154                                vc->Mail();
2155                               
2156                                if (setCounter)
2157                                        vc->mCounter = 0;
2158                                //viewCells.push_back(vc);
2159                        }
2160
2161                        viewCells.push_back(vc);
2162                }
2163
2164                if (setCounter)
2165                        ++ vc->mCounter;
2166        }
2167
2168        mCollectTimer.Exit();
2169
2170        return (int)objViewCells->size();
2171}
2172
2173
2174int BvHierarchy::CollectViewCellsFromRays(Intersectable *obj,
2175                                                                                  ViewCellContainer &viewCells,
2176                                                                                  const bool useMailBoxing,
2177                                                                                  const bool setCounter,
2178                                                                                  const bool onlyUnmailedRays)
2179{
2180        mCollectTimer.Entry();
2181        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();
2182
2183        int numRays = 0;
2184
2185        for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
2186        {
2187                VssRay *ray = (*rit);
2188
2189                if (onlyUnmailedRays && ray->Mailed())
2190                        continue;
2191               
2192                ++ numRays;
2193
2194                ViewCellContainer tmpViewCells;
2195                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
2196
2197                // matt: probably slow to allocate memory for view cells every time
2198                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
2199
2200                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
2201                {
2202                        ViewCell *vc = *vit;
2203
2204                        // store view cells
2205                        if (!useMailBoxing || !vc->Mailed())
2206                        {
2207                                if (useMailBoxing) // => view cell not mailed
2208                                {
2209                                        vc->Mail();
2210                                        if (setCounter)
2211                                                vc->mCounter = 0;
2212                                }
2213
2214                                viewCells.push_back(vc);
2215                        }
2216                       
2217                        if (setCounter)
2218                                ++ vc->mCounter;
2219                }
2220        }
2221
2222        mCollectTimer.Exit();
2223        return numRays;
2224}
2225
2226
2227int BvHierarchy::CountViewCellsFromRays(Intersectable *obj) //const
2228{
2229        int result = 0;
2230       
2231        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();
2232
2233        for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
2234        {
2235                VssRay *ray = (*rit);
2236                ViewCellContainer tmpViewCells;
2237       
2238                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
2239               
2240                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
2241                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
2242                {
2243                        ViewCell *vc = *vit;
2244
2245                        // store view cells
2246                        if (!vc->Mailed())
2247                        {
2248                                vc->Mail();
2249                                ++ result;
2250                        }
2251                }
2252        }
2253
2254        return result;
2255}
2256
2257#else
2258
2259int BvHierarchy::CountViewCells(Intersectable *obj) //const
2260{
2261        int result = 0;
2262       
2263        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();
2264
2265        for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
2266        {
2267                VssRay *ray = (*rit);
2268                ViewCellContainer tmpViewCells;
2269       
2270                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
2271               
2272                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
2273                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
2274                {
2275                        ViewCell *vc = *vit;
2276
2277                        // store view cells
2278                        if (!vc->Mailed())
2279                        {
2280                                vc->Mail();
2281                                ++ result;
2282                        }
2283                }
2284        }
2285
2286        return result;
2287}
2288
2289
2290int BvHierarchy::CollectViewCells(Intersectable *obj,
2291                                                                  ViewCellContainer &viewCells,
2292                                                                  const bool useMailBoxing,
2293                                                                  const bool setCounter,
2294                                                                  const bool onlyUnmailedRays)
2295{
2296        mCollectTimer.Entry();
2297        VssRayContainer::const_iterator rit, rit_end = obj->GetOrCreateRays()->end();
2298
2299        int numRays = 0;
2300
2301        for (rit = obj->GetOrCreateRays()->begin(); rit < rit_end; ++ rit)
2302        {
2303                VssRay *ray = (*rit);
2304
2305                if (onlyUnmailedRays && ray->Mailed())
2306                        continue;
2307               
2308                ++ numRays;
2309
2310                ViewCellContainer tmpViewCells;
2311                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
2312
2313                // matt: probably slow to allocate memory for view cells every time
2314                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
2315
2316                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
2317                {
2318                        ViewCell *vc = *vit;
2319
2320                        // store view cells
2321                        if (!useMailBoxing || !vc->Mailed())
2322                        {
2323                                if (useMailBoxing) // => view cell not mailed
2324                                {
2325                                        vc->Mail();
2326                                        if (setCounter)
2327                                                vc->mCounter = 0;
2328                                }
2329
2330                                viewCells.push_back(vc);
2331                        }
2332                       
2333                        if (setCounter)
2334                                ++ vc->mCounter;
2335                }
2336        }
2337
2338        mCollectTimer.Exit();
2339        return numRays;
2340}
2341#endif
2342
2343
2344int BvHierarchy::CountViewCells(const ObjectContainer &objects)// const
2345{
2346        int nViewCells = 0;
2347
2348        ViewCell::NewMail();
2349        ObjectContainer::const_iterator oit, oit_end = objects.end();
2350
2351        // loop through all object and collect view cell pvs of this node
2352        for (oit = objects.begin(); oit != oit_end; ++ oit)
2353        {
2354                nViewCells += CountViewCells(*oit);
2355        }
2356
2357        return nViewCells;
2358}
2359
2360
2361void BvHierarchy::CollectDirtyCandidates(BvhSubdivisionCandidate *sc,
2362                                                                                 vector<SubdivisionCandidate *> &dirtyList,
2363                                                                                 const bool onlyUnmailed)
2364{
2365        BvhTraversalData &tData = sc->mParentData;
2366        BvhLeaf *node = tData.mNode;
2367       
2368        ViewCellContainer viewCells;
2369        //ViewCell::NewMail();
2370        int numRays = CollectViewCells(*tData.mSampledObjects, viewCells, false, false);
2371
2372        if (0) cout << "collected " << (int)viewCells.size() << " dirty candidates" << endl;
2373       
2374        // split candidates handling
2375        // these view cells  are thrown into dirty list
2376        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
2377
2378        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
2379        {
2380        VspViewCell *vc = static_cast<VspViewCell *>(*vit);
2381                VspLeaf *leaf = vc->mLeaves[0];
2382       
2383                SubdivisionCandidate *candidate = leaf->GetSubdivisionCandidate();
2384               
2385                // is this leaf still a split candidate?
2386                if (candidate && (!onlyUnmailed || !candidate->Mailed()))
2387                {
2388                        candidate->Mail();
2389                        candidate->SetDirty(true);
2390                        dirtyList.push_back(candidate);
2391                }
2392        }
2393}
2394
2395
2396BvhNode *BvHierarchy::GetRoot() const
2397{
2398        return mRoot;
2399}
2400
2401
2402bool BvHierarchy::IsObjectInLeaf(BvhLeaf *leaf, Intersectable *object) const
2403{
2404        ObjectContainer::const_iterator oit =
2405                lower_bound(leaf->mObjects.begin(), leaf->mObjects.end(), object, ilt);
2406                               
2407        // objects sorted by id
2408        if ((oit != leaf->mObjects.end()) && ((*oit)->GetId() == object->GetId()))
2409        {
2410                return true;
2411        }
2412        else
2413        {
2414                return false;
2415        }
2416}
2417
2418#if 0
2419BvhLeaf *BvHierarchy::GetLeaf(Intersectable *object, BvhNode *node) const
2420{
2421        // hack: we use the simpler but faster version
2422        if (!object)
2423                return NULL;
2424
2425        return object->mBvhLeaf;
2426       
2427        ///////////////////////////////////////
2428        // start from root of tree
2429
2430        if (node == NULL)
2431                node = mRoot;
2432       
2433        vector<BvhLeaf *> leaves;
2434
2435        stack<BvhNode *> nodeStack;
2436        nodeStack.push(node);
2437 
2438        BvhLeaf *leaf = NULL;
2439 
2440        while (!nodeStack.empty()) 
2441        {
2442                BvhNode *node = nodeStack.top();
2443                nodeStack.pop();
2444       
2445                if (node->IsLeaf())
2446                {
2447                        leaf = static_cast<BvhLeaf *>(node);
2448
2449                        if (IsObjectInLeaf(leaf, object))
2450                        {
2451                                return leaf;
2452                        }
2453                }
2454                else   
2455                {       
2456                        // find point
2457                        BvhInterior *interior = static_cast<BvhInterior *>(node);
2458       
2459                        if (interior->GetBack()->GetBoundingBox().Includes(object->GetBox()))
2460                        {
2461                                nodeStack.push(interior->GetBack());
2462                        }
2463                       
2464                        // search both sides as we are using bounding volumes
2465                        if (interior->GetFront()->GetBoundingBox().Includes(object->GetBox()))
2466                        {
2467                                nodeStack.push(interior->GetFront());
2468                        }
2469                }
2470        }
2471 
2472        return leaf;
2473}
2474#endif
2475
2476bool BvHierarchy::Export(OUT_STREAM &stream)
2477{
2478        ExportNode(mRoot, stream);
2479
2480        return true;
2481}
2482
2483
2484void BvHierarchy::ExportObjects(BvhLeaf *leaf, OUT_STREAM &stream)
2485{
2486        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
2487
2488        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
2489        {
2490                stream << (*oit)->GetId() << " ";
2491        }
2492}
2493
2494
2495void BvHierarchy::ExportNode(BvhNode *node, OUT_STREAM &stream)
2496{
2497        if (node->IsLeaf())
2498        {
2499                BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
2500                const AxisAlignedBox3 box = leaf->GetBoundingBox();
2501                stream << "<Leaf id=\"" << node->GetId() << "\""
2502                           << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
2503                           << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\""
2504                           << " objects=\"";
2505               
2506                //-- export objects
2507                // tmp matt
2508                if (1) ExportObjects(leaf, stream);
2509               
2510                stream << "\" />" << endl;
2511        }
2512        else
2513        {       
2514                BvhInterior *interior = static_cast<BvhInterior *>(node);
2515                const AxisAlignedBox3 box = interior->GetBoundingBox();
2516
2517                stream << "<Interior id=\"" << node->GetId() << "\""
2518                           << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
2519                           << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z
2520                           << "\">" << endl;
2521
2522                ExportNode(interior->GetBack(), stream);
2523                ExportNode(interior->GetFront(), stream);
2524
2525                stream << "</Interior>" << endl;
2526        }
2527}
2528
2529
2530float BvHierarchy::EvalViewCellsVolume(const ObjectContainer &objects)// const
2531{
2532        float vol = 0;
2533
2534        ViewCellContainer viewCells;
2535       
2536        // we have to account for all view cells that can
2537        // be seen from the objects
2538        int numRays = CollectViewCells(objects, viewCells, false, false);
2539
2540        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
2541
2542        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
2543        {
2544                vol += (*vit)->GetVolume();
2545        }
2546
2547        return vol;
2548}
2549
2550
2551void BvHierarchy::Initialise(const ObjectContainer &objects)
2552{
2553        AxisAlignedBox3 box = EvalBoundingBox(objects);
2554
2555        ///////
2556        //-- create new root
2557
2558        BvhLeaf *bvhleaf = new BvhLeaf(box, NULL, (int)objects.size());
2559        bvhleaf->mObjects = objects;
2560        mRoot = bvhleaf;
2561
2562        // compute bounding box from objects
2563        mBoundingBox = mRoot->GetBoundingBox();
2564
2565        // associate root with current objects
2566        AssociateObjectsWithLeaf(bvhleaf);
2567}
2568
2569
2570void BvHierarchy::StoreSampledObjects(ObjectContainer &sampledObjects, const ObjectContainer &objects)
2571{
2572        ObjectContainer::const_iterator oit, oit_end = objects.end();
2573
2574        for (oit = objects.begin(); oit != objects.end(); ++ oit)
2575        {
2576                Intersectable *obj = *oit;
2577
2578                if (!obj->GetOrCreateRays()->empty())
2579        {
2580                        sampledObjects.push_back(obj);
2581                }
2582        }
2583        }
2584
2585
2586void BvHierarchy::PrepareConstruction(SplitQueue &tQueue,
2587                                                                          const VssRayContainer &sampleRays,
2588                                                                          const ObjectContainer &objects)
2589{
2590        ///////////////////////////////////////
2591        //-- we assume that we have objects sorted by their id =>
2592        //-- we don't have to sort them here and an binary search
2593        //-- for identifying if a object is in a leaf.
2594       
2595        mBvhStats.Reset();
2596        mBvhStats.Start();
2597        mBvhStats.nodes = 1;
2598               
2599        // store pointer to this tree
2600        BvhSubdivisionCandidate::sBvHierarchy = this;
2601       
2602        // root and bounding box was already constructed
2603        BvhLeaf *bvhLeaf = static_cast<BvhLeaf *>(mRoot);
2604       
2605        // only rays intersecting objects in node are interesting
2606        const int nRays = AssociateObjectsWithRays(sampleRays);
2607        //cout << "using " << nRays << " of " << (int)sampleRays.size() << " rays" << endl;
2608       
2609        ObjectContainer *sampledObjects = new ObjectContainer();
2610        StoreSampledObjects(*sampledObjects, objects);
2611
2612#if STORE_VIEWCELLS_WITH_BVH
2613        AssociateViewCellsWithObjects(*sampledObjects);
2614#endif
2615
2616        // probability that volume is "seen" from the view cells
2617        const float prop = EvalViewCellsVolume(*sampledObjects) / GetViewSpaceVolume();
2618
2619        // create bvh traversal data
2620        BvhTraversalData oData(bvhLeaf, 0, prop, nRays);
2621               
2622        // create sorted object lists for the first data
2623        if (mUseGlobalSorting)
2624        {
2625                AssignInitialSortedObjectList(oData, objects);
2626        }
2627       
2628        oData.mSampledObjects = sampledObjects;
2629       
2630        ///////////////////
2631        //-- add first candidate for object space partition     
2632
2633        mTotalCost = EvalRenderCost(objects);
2634        mPvsEntries = CountViewCells(*sampledObjects);
2635
2636        oData.mCorrectedPvs = oData.mPvs = (float)mPvsEntries;
2637        oData.mCorrectedVolume = oData.mVolume = prop;
2638       
2639        BvhSubdivisionCandidate *oSubdivisionCandidate =
2640                new BvhSubdivisionCandidate(oData);
2641
2642        bvhLeaf->SetSubdivisionCandidate(oSubdivisionCandidate);
2643
2644#if STORE_VIEWCELLS_WITH_BVH
2645        ReleaseViewCells(*sampledObjects);
2646#endif
2647
2648        if (mApplyInitialPartition)
2649        {
2650                vector<SubdivisionCandidate *> candidateContainer;
2651
2652                mIsInitialSubdivision = true;
2653               
2654                // evaluate priority
2655                EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
2656                PrintSubdivisionStats(*oSubdivisionCandidate);
2657
2658                ApplyInitialSubdivision(oSubdivisionCandidate, candidateContainer);             
2659
2660                mIsInitialSubdivision = false;
2661
2662                vector<SubdivisionCandidate *>::const_iterator cit, cit_end = candidateContainer.end();
2663
2664                for (cit = candidateContainer.begin(); cit != cit_end; ++ cit)
2665                {
2666                        BvhSubdivisionCandidate *sCandidate = static_cast<BvhSubdivisionCandidate *>(*cit);
2667                       
2668                        // reevaluate priority
2669                        EvalSubdivisionCandidate(*sCandidate, true, true);
2670                        tQueue.Push(sCandidate);
2671                }
2672
2673                cout << "size of initial bv subdivision: " << GetStatistics().Leaves() << endl;
2674        }
2675        else
2676        {       
2677                // evaluate priority
2678                EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
2679                PrintSubdivisionStats(*oSubdivisionCandidate);
2680
2681                tQueue.Push(oSubdivisionCandidate);
2682                cout << "size of initial bv subdivision: " << GetStatistics().Leaves() << endl;
2683        }
2684}
2685
2686
2687void BvHierarchy::AssignInitialSortedObjectList(BvhTraversalData &tData,
2688                                                                                                const ObjectContainer &objects)
2689{
2690        const bool doSort = true;
2691
2692        // we sort the objects as a preprocess so they don't have
2693        // to be sorted for each split
2694        for (int i = 0; i < 3; ++ i)
2695        {
2696                SortableEntryContainer *sortedObjects = new SortableEntryContainer();
2697
2698                CreateLocalSubdivisionCandidates(objects,
2699                                                                             &sortedObjects,
2700                                                                                 doSort,
2701                                                                                 i);
2702               
2703                // copy list into traversal data list
2704                tData.mSortedObjects[i] = new ObjectContainer();
2705                tData.mSortedObjects[i]->reserve((int)objects.size());
2706
2707                SortableEntryContainer::const_iterator oit, oit_end = sortedObjects->end();
2708
2709                for (oit = sortedObjects->begin(); oit != oit_end; ++ oit)
2710                {
2711                        tData.mSortedObjects[i]->push_back((*oit).mObject);
2712                }
2713
2714                delete sortedObjects;
2715        }
2716
2717        // next sorted list: by size (for initial criteria)
2718        tData.mSortedObjects[3] = new ObjectContainer();
2719        tData.mSortedObjects[3]->reserve((int)objects.size());
2720
2721        *(tData.mSortedObjects[3]) = objects;
2722       
2723        stable_sort(tData.mSortedObjects[3]->begin(), tData.mSortedObjects[3]->end(), smallerSize);
2724}
2725
2726
2727void BvHierarchy::AssignSortedObjects(const BvhSubdivisionCandidate &sc,
2728                                                                          BvhTraversalData &frontData,
2729                                                                          BvhTraversalData &backData)
2730{
2731        Intersectable::NewMail();
2732
2733        // we sorted the objects as a preprocess so they don't have
2734        // to be sorted for each split
2735        ObjectContainer::const_iterator fit, fit_end = sc.mFrontObjects.end();
2736
2737        for (fit = sc.mFrontObjects.begin(); fit != fit_end; ++ fit)
2738        {
2739                (*fit)->Mail();
2740        }
2741
2742        for (int i = 0; i < 4; ++ i)
2743        {
2744                frontData.mSortedObjects[i] = new ObjectContainer();
2745                backData.mSortedObjects[i] = new ObjectContainer();
2746
2747                frontData.mSortedObjects[i]->reserve((int)sc.mFrontObjects.size());
2748                backData.mSortedObjects[i]->reserve((int)sc.mBackObjects.size());
2749
2750                ObjectContainer::const_iterator oit, oit_end = sc.mParentData.mSortedObjects[i]->end();
2751
2752                // all the front objects are mailed => assign the sorted object lists
2753                for (oit = sc.mParentData.mSortedObjects[i]->begin(); oit != oit_end; ++ oit)
2754                {
2755                        if ((*oit)->Mailed())
2756                        {
2757                                frontData.mSortedObjects[i]->push_back(*oit);
2758                        }
2759                        else
2760                        {
2761                                backData.mSortedObjects[i]->push_back(*oit);
2762                        }
2763                }
2764        }
2765}
2766
2767
2768void BvHierarchy::Reset(SplitQueue &tQueue,
2769                                                const VssRayContainer &sampleRays,
2770                                                const ObjectContainer &objects)
2771{
2772
2773        // reset stats
2774        mBvhStats.Reset();
2775        mBvhStats.Start();
2776        mBvhStats.nodes = 1;
2777
2778        // reset root
2779        DEL_PTR(mRoot);
2780       
2781        BvhLeaf *bvhleaf = new BvhLeaf(mBoundingBox, NULL, (int)objects.size());
2782        bvhleaf->mObjects = objects;
2783        mRoot = bvhleaf;
2784       
2785        ObjectContainer *sampledObjects = new ObjectContainer();
2786        StoreSampledObjects(*sampledObjects, objects);
2787
2788#if STORE_VIEWCELLS_WITH_BVH
2789        AssociateViewCellsWithObjects(*sampledObjects);
2790#endif
2791
2792        //mTermMinProbability *= mVspTree->GetBoundingBox().GetVolume();
2793        // probability that volume is "seen" from the view cells
2794        const float viewSpaceVol = mViewCellsManager->GetViewSpaceBox().GetVolume();
2795        const float prop = EvalViewCellsVolume(*sampledObjects);
2796
2797        const int nRays = CountRays(*sampledObjects);
2798        BvhLeaf *bvhLeaf = static_cast<BvhLeaf *>(mRoot);
2799
2800        // create bvh traversal data
2801        BvhTraversalData oData(bvhLeaf, 0, prop, nRays);
2802
2803        oData.mSampledObjects = sampledObjects;
2804
2805        if (mUseGlobalSorting)
2806                AssignInitialSortedObjectList(oData, objects);
2807       
2808#if STORE_VIEWCELLS_WITH_BVH
2809        ReleaseViewCells(*sampledObjects);
2810#endif
2811        ///////////////////
2812        //-- add first candidate for object space partition     
2813
2814        BvhSubdivisionCandidate *oSubdivisionCandidate =
2815                new BvhSubdivisionCandidate(oData);
2816
2817        EvalSubdivisionCandidate(*oSubdivisionCandidate, true, true);
2818        bvhLeaf->SetSubdivisionCandidate(oSubdivisionCandidate);
2819
2820        mTotalCost = (float)objects.size() * prop;
2821
2822        PrintSubdivisionStats(*oSubdivisionCandidate);
2823
2824        tQueue.Push(oSubdivisionCandidate);
2825}
2826
2827
2828void BvhStatistics::Print(ostream &app) const
2829{
2830        app << "=========== BvHierarchy statistics ===============\n";
2831
2832        app << setprecision(4);
2833
2834        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
2835
2836        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
2837
2838        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
2839
2840        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
2841
2842        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits << endl;
2843
2844        app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
2845                << maxCostNodes * 100 / (double)Leaves() << endl;
2846
2847        app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
2848                << minProbabilityNodes * 100 / (double)Leaves() << endl;
2849
2850
2851        //////////////////////////////////////////////////
2852       
2853        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
2854                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
2855       
2856        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
2857
2858        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
2859
2860        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
2861
2862       
2863        ////////////////////////////////////////////////////////
2864       
2865        app << "#N_PMINOBJECTSLEAVES  ( Percentage of leaves with mininum objects )\n"
2866                << minObjectsNodes * 100 / (double)Leaves() << endl;
2867
2868        app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" << maxObjectRefs << "\n";
2869
2870        app << "#N_MINOBJECTREFS  ( Min number of object refs / leaf )\n" << minObjectRefs << "\n";
2871
2872        app << "#N_EMPTYLEAFS ( Empty leafs )\n" << emptyNodes << "\n";
2873       
2874        app << "#N_PAVGOBJECTSLEAVES  ( average object refs / leaf)\n" << AvgObjectRefs() << endl;
2875
2876
2877        ////////////////////////////////////////////////////////
2878       
2879        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with mininum rays )\n"
2880                << minRaysNodes * 100 / (double)Leaves() << endl;
2881
2882        app << "#N_MAXRAYREFS  ( Max number of ray refs / leaf )\n" << maxRayRefs << "\n";
2883
2884        app << "#N_MINRAYREFS  ( Min number of ray refs / leaf )\n" << minRayRefs << "\n";
2885       
2886        app << "#N_PAVGRAYLEAVES  ( average ray refs / leaf )\n" << AvgRayRefs() << endl;
2887       
2888        app << "#N_PAVGRAYCONTRIBLEAVES  ( Average ray contribution)\n" <<
2889                rayRefs / (double)objectRefs << endl;
2890
2891        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
2892                maxRayContriNodes * 100 / (double)Leaves() << endl;
2893
2894        app << "#N_PGLOBALCOSTMISSES ( Global cost misses )\n" << mGlobalCostMisses << endl;
2895
2896        app << "========== END OF BvHierarchy statistics ==========\n";
2897}
2898
2899
2900// TODO: return memory usage in MB
2901float BvHierarchy::GetMemUsage() const
2902{
2903        return (float)(sizeof(BvHierarchy)
2904                                   + mBvhStats.Leaves() * sizeof(BvhLeaf)
2905                                   + mBvhStats.Interior() * sizeof(BvhInterior)
2906                                   ) / float(1024 * 1024);
2907}
2908
2909
2910void BvHierarchy::SetActive(BvhNode *node) const
2911{
2912        vector<BvhLeaf *> leaves;
2913
2914        // sets the pointers to the currently active view cells
2915        CollectLeaves(node, leaves);
2916        vector<BvhLeaf *>::const_iterator lit, lit_end = leaves.end();
2917
2918        for (lit = leaves.begin(); lit != lit_end; ++ lit)
2919        {
2920                (*lit)->SetActiveNode(node);
2921        }
2922}
2923
2924
2925BvhNode *BvHierarchy::SubdivideAndCopy(SplitQueue &tQueue,
2926                                                                           SubdivisionCandidate *splitCandidate)
2927{
2928        BvhSubdivisionCandidate *sc =
2929                static_cast<BvhSubdivisionCandidate *>(splitCandidate);
2930        BvhTraversalData &tData = sc->mParentData;
2931
2932        BvhNode *currentNode = tData.mNode;
2933        BvhNode *oldNode = (BvhNode *)splitCandidate->mEvaluationHack;
2934
2935        if (!oldNode->IsLeaf())
2936        {       
2937                //////////////
2938                //-- continue subdivision
2939
2940                BvhTraversalData tFrontData;
2941                BvhTraversalData tBackData;
2942                       
2943                BvhInterior *oldInterior = static_cast<BvhInterior *>(oldNode);
2944               
2945                sc->mFrontObjects.clear();
2946                sc->mBackObjects.clear();
2947
2948                oldInterior->GetFront()->CollectObjects(sc->mFrontObjects);
2949                oldInterior->GetBack()->CollectObjects(sc->mBackObjects);
2950               
2951                // evaluate the changes in render cost and pvs entries
2952                EvalSubdivisionCandidate(*sc, false, true);
2953
2954                // create new interior node and two leaf node
2955                currentNode = SubdivideNode(*sc, tFrontData, tBackData);
2956       
2957                //oldNode->mRenderCostDecr += sc->GetRenderCostDecrease();
2958                //oldNode->mPvsEntriesIncr += sc->GetPvsEntriesIncr();
2959               
2960                //oldNode->mRenderCostDecr = sc->GetRenderCostDecrease();
2961                //oldNode->mPvsEntriesIncr = sc->GetPvsEntriesIncr();
2962               
2963                ///////////////////////////
2964                //-- push the new split candidates on the queue
2965               
2966                BvhSubdivisionCandidate *frontCandidate = new BvhSubdivisionCandidate(tFrontData);
2967                BvhSubdivisionCandidate *backCandidate = new BvhSubdivisionCandidate(tBackData);
2968
2969                frontCandidate->SetPriority((float)-oldInterior->GetFront()->GetTimeStamp());
2970                backCandidate->SetPriority((float)-oldInterior->GetBack()->GetTimeStamp());
2971
2972                frontCandidate->mEvaluationHack = oldInterior->GetFront();
2973                backCandidate->mEvaluationHack = oldInterior->GetBack();
2974
2975                // cross reference
2976                tFrontData.mNode->SetSubdivisionCandidate(frontCandidate);
2977                tBackData.mNode->SetSubdivisionCandidate(backCandidate);
2978
2979                //cout << "f: " << frontCandidate->GetPriority() << " b: " << backCandidate->GetPriority() << endl;
2980                tQueue.Push(frontCandidate);
2981                tQueue.Push(backCandidate);
2982        }
2983
2984        /////////////////////////////////
2985        //-- node is a leaf => terminate traversal
2986
2987        if (currentNode->IsLeaf())
2988        {
2989                // this leaf is no candidate for splitting anymore
2990                // => detach subdivision candidate
2991                tData.mNode->SetSubdivisionCandidate(NULL);
2992                // detach node so we don't delete it with the traversal data
2993                tData.mNode = NULL;
2994        }
2995       
2996        return currentNode;
2997}
2998
2999
3000void BvHierarchy::CollectObjects(const AxisAlignedBox3 &box,
3001                                                                 ObjectContainer &objects)
3002{
3003  stack<BvhNode *> nodeStack;
3004 
3005  nodeStack.push(mRoot);
3006
3007  while (!nodeStack.empty()) {
3008        BvhNode *node = nodeStack.top();
3009       
3010        nodeStack.pop();
3011        if (node->IsLeaf()) {
3012          BvhLeaf *leaf = (BvhLeaf *)node;
3013          if (Overlap(box, leaf->GetBoundingBox())) {
3014                Intersectable *object = leaf;
3015                if (!object->Mailed()) {
3016                  object->Mail();
3017                  objects.push_back(object);
3018                }
3019          }
3020        }
3021        else
3022          {
3023                BvhInterior *interior = (BvhInterior *)node;
3024                if (Overlap(box, interior->GetBoundingBox())) {
3025                  bool pushed = false;
3026                  if (!interior->GetFront()->Mailed()) {
3027                        nodeStack.push(interior->GetFront());
3028                        pushed = true;
3029                  }
3030                  if (!interior->GetBack()->Mailed()) {
3031                        nodeStack.push(interior->GetBack());
3032                        pushed = true;
3033                  }
3034                  // avoid traversal of this node in the next query
3035                  if (!pushed)
3036                        interior->Mail();
3037                }
3038          }
3039  }
3040}
3041
3042
3043void BvHierarchy::CreateUniqueObjectIds()
3044{
3045        stack<BvhNode *> nodeStack;
3046        nodeStack.push(mRoot);
3047
3048        int currentId = 0;
3049        while (!nodeStack.empty())
3050        {
3051                BvhNode *node = nodeStack.top();
3052                nodeStack.pop();
3053
3054                node->SetId(currentId ++);
3055
3056                if (!node->IsLeaf())
3057                {
3058                        BvhInterior *interior = (BvhInterior *)node;
3059
3060                        nodeStack.push(interior->GetFront());
3061                        nodeStack.push(interior->GetBack());
3062                }
3063        }
3064}
3065
3066
3067void BvHierarchy::ApplyInitialSubdivision(SubdivisionCandidate *firstCandidate,
3068                                                                                  vector<SubdivisionCandidate *> &candidateContainer)
3069{
3070        SplitQueue tempQueue;
3071        tempQueue.Push(firstCandidate);
3072
3073        while (!tempQueue.Empty())
3074        {
3075                SubdivisionCandidate *candidate = tempQueue.Top();
3076                tempQueue.Pop();
3077
3078                BvhSubdivisionCandidate *bsc =
3079                        static_cast<BvhSubdivisionCandidate *>(candidate);
3080
3081                if (!InitialTerminationCriteriaMet(bsc->mParentData))
3082                {
3083                        const bool globalCriteriaMet = GlobalTerminationCriteriaMet(bsc->mParentData);
3084               
3085                        SubdivisionCandidateContainer dirtyList;
3086                        BvhNode *node = Subdivide(tempQueue, bsc, globalCriteriaMet, dirtyList);
3087
3088                        // not needed anymore
3089                        delete bsc;
3090                }
3091                else
3092                {
3093                        // initial preprocessing  finished for this candidate
3094                        // add to candidate container
3095                        candidateContainer.push_back(bsc);
3096                }
3097        }
3098}
3099
3100
3101void BvHierarchy::ApplyInitialSplit(const BvhTraversalData &tData,
3102                                                                        ObjectContainer &frontObjects,
3103                                                                        ObjectContainer &backObjects)
3104{
3105        ObjectContainer *objects = tData.mSortedObjects[3];
3106
3107        ObjectContainer::const_iterator oit, oit_end = objects->end();
3108   
3109        float maxAreaDiff = -1.0f;
3110
3111        ObjectContainer::const_iterator backObjectsStart = objects->begin();
3112
3113        for (oit = objects->begin(); oit != (objects->end() - 1); ++ oit)
3114        {
3115                Intersectable *objS = *oit;
3116                Intersectable *objL = *(oit + 1);
3117               
3118                const float areaDiff =
3119                                objL->GetBox().SurfaceArea() - objS->GetBox().SurfaceArea();
3120
3121                if (areaDiff > maxAreaDiff)
3122                {
3123                        maxAreaDiff = areaDiff;
3124                        backObjectsStart = oit + 1;
3125                }
3126        }
3127
3128        // belongs to back bv
3129        for (oit = objects->begin(); oit != backObjectsStart; ++ oit)
3130        {
3131                frontObjects.push_back(*oit);
3132        }
3133
3134        // belongs to front bv
3135        for (oit = backObjectsStart; oit != oit_end; ++ oit)
3136        {
3137                backObjects.push_back(*oit);
3138        }
3139       
3140        cout << "front: " << (int)frontObjects.size() << " back: " << (int)backObjects.size() << " "
3141                 << backObjects.front()->GetBox().SurfaceArea() - frontObjects.back()->GetBox().SurfaceArea() << endl;
3142}
3143
3144
3145inline static float AreaRatio(Intersectable *smallObj, Intersectable *largeObj)
3146{
3147        const float areaSmall = smallObj->GetBox().SurfaceArea();
3148        const float areaLarge = largeObj->GetBox().SurfaceArea();
3149
3150        return areaSmall / (areaLarge - areaSmall + Limits::Small);
3151}
3152
3153
3154bool BvHierarchy::InitialTerminationCriteriaMet(const BvhTraversalData &tData) const
3155{
3156        const bool terminationCriteriaMet =
3157                        (0
3158                    || ((int)tData.mNode->mObjects.size() < mInitialMinObjects)
3159                        || (tData.mNode->mObjects.back()->GetBox().SurfaceArea() < mInitialMinArea)
3160                        || (AreaRatio(tData.mNode->mObjects.front(), tData.mNode->mObjects.back()) > mInitialMaxAreaRatio)
3161                        );
3162
3163        cout << "criteria met: "<< terminationCriteriaMet << "\n"
3164                 << "size: " << (int)tData.mNode->mObjects.size() << " max: " << mInitialMinObjects << endl
3165                 << "ratio: " << AreaRatio(tData.mNode->mObjects.front(), tData.mNode->mObjects.back()) << " max: " << mInitialMaxAreaRatio << endl
3166                 << "area: " << tData.mNode->mObjects.back()->GetBox().SurfaceArea() << " max: " << mInitialMinArea << endl << endl;
3167
3168        return terminationCriteriaMet;
3169}
3170
3171
3172// HACK
3173float BvHierarchy::GetRenderCostIncrementially(BvhNode *node) const
3174{
3175        if (node->mRenderCost < 0) // render cost must be evaluated ...
3176        {
3177                if (node->IsLeaf())
3178                {
3179                        BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
3180                        node->mRenderCost = EvalAbsCost(leaf->mObjects);               
3181                        //if (leaf->mObjects.size() < 100) cout <<"p";
3182                }
3183                else
3184                {
3185                        BvhInterior *interior = static_cast<BvhInterior *>(node);
3186               
3187                        node->mRenderCost = GetRenderCostIncrementially(interior->GetFront()) +
3188                                                                GetRenderCostIncrementially(interior->GetBack());
3189                }
3190        }
3191
3192        return node->mRenderCost;
3193}
3194
3195
3196void BvHierarchy::Compress()
3197{
3198}
3199
3200
3201void BvHierarchy::SetUniqueNodeIds()
3202{
3203        // export bounding boxes
3204        vector<BvhNode *> nodes;
3205
3206        // hack: should also expect interior nodes
3207        CollectNodes(mRoot, nodes);
3208
3209        vector<BvhNode *>::const_iterator oit, oit_end = nodes.end();
3210
3211        int id = 0;
3212
3213        for (oit = nodes.begin(); oit != oit_end; ++ oit, ++ id)
3214        {
3215                (*oit)->SetId(id);
3216        }
3217}
3218
3219
3220}
Note: See TracBrowser for help on using the repository browser.