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

Revision 2588, 83.2 KB checked in by bittner, 17 years ago (diff)

sceneBox bugfix

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