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

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