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

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