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

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