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

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