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

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