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

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