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

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