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

Revision 1414, 48.8 KB checked in by mattausch, 18 years ago (diff)

fixed kd tree loading / exporting

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
297        mBvhStats.nodes += 2; // we have two new leaves
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        // matt: TODO
525        return (0
526                || (mBvhStats.Leaves() >= mTermMaxLeaves)
527                || (mGlobalCostMisses >= mTermGlobalCostMissTolerance)
528                //|| mOutOfMemory
529                );
530}
531
532
533void BvHierarchy::EvaluateLeafStats(const BvhTraversalData &data)
534{
535        // the node became a leaf -> evaluate stats for leafs
536        BvhLeaf *leaf = data.mNode;
537       
538        ++ mCreatedLeaves;
539
540       
541        if (data.mProbability <= mTermMinProbability)
542        {
543                ++ mBvhStats.minProbabilityNodes;
544        }
545
546        ////////////////////////////////////////////
547        // depth related stuff
548
549        if (data.mDepth < mBvhStats.minDepth)
550        {
551                mBvhStats.minDepth = data.mDepth;
552        }
553
554        if (data.mDepth >= mTermMaxDepth)
555        {
556        ++ mBvhStats.maxDepthNodes;
557        }
558
559        // accumulate depth to compute average depth
560        mBvhStats.accumDepth += data.mDepth;
561
562
563        ////////////////////////////////////////////
564        // objects related stuff
565
566        // note: this number should always accumulate to the total number of objects
567        mBvhStats.objectRefs += (int)leaf->mObjects.size();
568
569        if ((int)leaf->mObjects.size() <= mTermMinObjects)
570        {
571             ++ mBvhStats.minObjectsNodes;
572        }
573
574        if (leaf->mObjects.empty())
575        {
576                ++ mBvhStats.emptyNodes;
577        }
578
579        if ((int)leaf->mObjects.size() > mBvhStats.maxObjectRefs)
580        {
581                mBvhStats.maxObjectRefs = (int)leaf->mObjects.size();
582        }
583
584        if ((int)leaf->mObjects.size() < mBvhStats.minObjectRefs)
585        {
586                mBvhStats.minObjectRefs = (int)leaf->mObjects.size();
587        }
588
589        ////////////////////////////////////////////
590        // ray related stuff
591
592        // note: this number should always accumulate to the total number of rays
593        mBvhStats.rayRefs += data.mNumRays;
594       
595        if (data.mNumRays <= mTermMinRays)
596        {
597             ++ mBvhStats.minRaysNodes;
598        }
599
600        if (data.mNumRays > mBvhStats.maxRayRefs)
601        {
602                mBvhStats.maxRayRefs = data.mNumRays;
603        }
604
605        if (data.mNumRays < mBvhStats.minRayRefs)
606        {
607                mBvhStats.minRayRefs = data.mNumRays;
608        }
609
610        cout << "depth: " << data.mDepth << " objects: " << (int)leaf->mObjects.size()
611                 << " rays: " << data.mNumRays << " rays / objects "
612                 << (float)data.mNumRays / (float)leaf->mObjects.size() << endl;
613}
614
615
616#if 0
617
618/// compute object boundaries using spatial mid split
619float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
620                                                                                        const int axis,
621                                                                                        ObjectContainer &objectsFront,
622                                                                                        ObjectContainer &objectsBack)
623{
624        const float maxBox = tData.mBoundingBox.Max(axis);
625        const float minBox = tData.mBoundingBox.Min(axis);
626
627        float midPoint = (maxBox + minBox) * 0.5f;
628
629        ObjectContainer::const_iterator oit, oit_end = tData.mNode->mObjects.end();
630       
631        for (oit = tData.mNode->mObjects.begin(); oit != oit_end; ++ oit)
632        {
633                Intersectable *obj = *oit;
634                const AxisAlignedBox3 box = obj->GetBox();
635
636                const float objMid = (box.Max(axis) + box.Min(axis)) * 0.5f;
637
638                // object mailed => belongs to back objects
639                if (objMid < midPoint)
640                {
641                        objectsBack.push_back(obj);
642                }
643                else
644                {
645                        objectsFront.push_back(obj);
646                }
647        }
648
649        const float oldRenderCost = EvalRenderCost(tData.mNode->mObjects);
650        const float newRenderCost =
651                EvalRenderCost(objectsFront) * EvalRenderCost(objectsBack);
652
653        const float ratio = newRenderCost / oldRenderCost;
654        return ratio;
655}
656
657#else
658
659/// compute object partition by getting balanced objects on the left and right side
660float BvHierarchy::EvalLocalObjectPartition(const BvhTraversalData &tData,
661                                                                                        const int axis,
662                                                                                        ObjectContainer &objectsFront,
663                                                                                        ObjectContainer &objectsBack)
664{
665        PrepareLocalSubdivisionCandidates(tData, axis);
666       
667        SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();
668
669        int i = 0;
670        const int border = (int)tData.mNode->mObjects.size() / 2;
671
672    for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ i)
673        {
674                Intersectable *obj = (*cit).mObject;
675
676                // object mailed => belongs to back objects
677                if (i < border)
678                {
679                        objectsBack.push_back(obj);
680                }
681                else
682                {
683                        objectsFront.push_back(obj);
684                }
685        }
686
687        const float oldRenderCost = EvalRenderCost(tData.mNode->mObjects);
688        const float newRenderCost = EvalRenderCost(objectsFront) + EvalRenderCost(objectsBack);
689
690        const float ratio = newRenderCost / oldRenderCost;
691        return ratio;
692}
693#endif
694
695
696float BvHierarchy::EvalSah(const BvhTraversalData &tData,
697                                                   const int axis,
698                                                   ObjectContainer &objectsFront,
699                                                   ObjectContainer &objectsBack)
700{
701        // go through the lists, count the number of objects left and right
702        // and evaluate the following cost funcion:
703        // C = ct_div_ci  + (ol + or)/queries
704        PrepareLocalSubdivisionCandidates(tData, axis);
705
706        int objectsLeft = 0, objectsRight = (int)tData.mNode->mObjects.size();
707 
708        AxisAlignedBox3 box = tData.mNode->GetBoundingBox();
709
710        float minBox = box.Min(axis);
711        float maxBox = box.Max(axis);
712        float boxArea = box.SurfaceArea();
713
714        float minSum = 1e20f;
715 
716        float minBorder = maxBox;
717        float maxBorder = minBox;
718        float areaLeft = 0, areaRight = 0;
719
720        SortableEntryContainer::const_iterator currentPos =
721                mSubdivisionCandidates->begin();
722
723       
724        // we keep track of both borders of the bounding boxes =>
725        // store the events in descending order
726        vector<float> bordersRight;
727        bordersRight.resize(mSubdivisionCandidates->size());
728
729        SortableEntryContainer::reverse_iterator rcit =
730                mSubdivisionCandidates->rbegin(), rcit_end = mSubdivisionCandidates->rend();
731       
732        vector<float>::reverse_iterator rbit = bordersRight.rbegin();
733       
734        for (; rcit != rcit_end; ++ rcit, ++ rbit)
735        {
736                Intersectable *obj = (*rcit).mObject;
737                const AxisAlignedBox3 box = obj->GetBox();
738
739                if (box.Min(axis) < minBorder)
740                {
741                        minBorder = box.Min(axis);
742                }
743
744                (*rbit) = minBorder;
745        }
746
747        vector<float>::const_iterator bit = bordersRight.begin();
748        SortableEntryContainer::const_iterator cit, cit_end = mSubdivisionCandidates->end();
749
750        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit, ++ bit)
751        {
752                Intersectable *obj = (*cit).mObject;
753
754                ++ objectsLeft;
755                -- objectsRight;
756
757                AxisAlignedBox3 lbox = box;
758                AxisAlignedBox3 rbox = box;
759       
760                const AxisAlignedBox3 obox = obj->GetBox();
761
762                // the borders of the bounding boxes have changed
763                if (obox.Max(axis) > maxBorder)
764                {
765                        maxBorder = obox.Max(axis);
766                }
767
768                minBorder = (*bit);
769       
770        lbox.SetMax(axis, maxBorder);
771                rbox.SetMin(axis, minBorder);
772       
773                const float al = lbox.SurfaceArea();
774                const float ar = rbox.SurfaceArea();
775
776                const float sum = objectsLeft * al + objectsRight * ar;
777     
778                /*cout << "pos=" << (*cit).mPos << "\t q=(" << objectsLeft << "," << objectsRight <<")\t r=("
779                         << lbox.SurfaceArea() << "," << rbox.SurfaceArea() << ")" << endl;
780                cout << "minborder: " << minBorder << " maxborder: " << maxBorder << endl;
781            cout << "cost= " << sum << endl;
782*/
783                if (sum < minSum)
784                {
785                        minSum = sum;
786                        areaLeft = al;
787                        areaRight = ar;
788                        // objects belong to left side now
789                        for (; currentPos != (cit + 1); ++ currentPos);
790                }
791        }
792
793
794        ////////////////////////////////////////////
795        //-- assign object to front and back volume
796
797        // belongs to back bv
798        for (cit = mSubdivisionCandidates->begin(); cit != currentPos; ++ cit)
799                objectsBack.push_back((*cit).mObject);
800       
801        // belongs to front bv
802        for (cit = currentPos; cit != cit_end; ++ cit)
803                objectsFront.push_back((*cit).mObject);
804
805        float oldCost = (float)tData.mNode->mObjects.size();
806        float newCost = minSum / boxArea;
807        float ratio = newCost / oldCost;
808 
809#ifdef _DEBUG
810        cout << "\n\nobjects=(" << (int)objectsBack.size() << "," << (int)objectsFront.size() << " of "
811                 << (int)tData.mNode->mObjects.size() << ")\t area=("
812                 << areaLeft << "," << areaRight << ")" << endl;
813        cout << "cost= " << minSum << endl;
814#endif
815  return ratio;
816}
817
818
819static bool PrepareOutput(const int axis,
820                                                  const int leaves,
821                                                  ofstream &sumStats,
822                                                  ofstream &vollStats,
823                                                  ofstream &volrStats)
824{
825        if ((axis == 0) && (leaves > 0) && (leaves < 90))
826        {
827                char str[64];   
828                sprintf(str, "tmp/bvh_heur_sum-%04d.log", leaves);
829                sumStats.open(str);
830                sprintf(str, "tmp/bvh_heur_voll-%04d.log", leaves);
831                vollStats.open(str);
832                sprintf(str, "tmp/bvh_heur_volr-%04d.log", leaves);
833                volrStats.open(str);
834        }
835
836        return sumStats.is_open() && vollStats.is_open() && volrStats.is_open();
837}
838
839
840static void PrintHeuristics(const int objectsRight,
841                                                        const float sum,
842                                                        const float volLeft,
843                                                        const float volRight,
844                                                        const float viewSpaceVol,
845                                                        ofstream &sumStats,
846                                                        ofstream &vollStats,
847                                                        ofstream &volrStats)
848{
849        sumStats
850                << "#Position\n" << objectsRight << endl
851                << "#Sum\n" << sum / viewSpaceVol << endl
852                << "#Vol\n" << (volLeft +  volRight) / viewSpaceVol << endl;
853
854        vollStats
855                << "#Position\n" << objectsRight << endl
856                << "#Vol\n" << volLeft / viewSpaceVol << endl;
857
858        volrStats
859                << "#Position\n" << objectsRight << endl
860                << "#Vol\n" << volRight / viewSpaceVol << endl;
861}
862
863
864float BvHierarchy::EvalLocalCostHeuristics(const BvhTraversalData &tData,
865                                                                                   const int axis,
866                                                                                   ObjectContainer &objectsFront,
867                                                                                   ObjectContainer &objectsBack)
868{
869        ////////////////////////////////////////////////////////////////
870        // go through the lists, count the number of objects left and right
871        // and evaluate the cost funcion
872
873        // prepare the heuristics by setting mailboxes and counters.
874        const float totalVol = PrepareHeuristics(tData, axis);
875       
876        // local helper variables
877        float volLeft = 0;
878        float volRight = totalVol;
879        int nObjectsLeft = 0;
880        const int nTotalObjects = (int)tData.mNode->mObjects.size();
881        const float viewSpaceVol = mHierarchyManager->GetViewSpaceBox().GetVolume();
882
883        SortableEntryContainer::const_iterator backObjectsStart = mSubdivisionCandidates->begin();
884
885        /////////////////////////////////
886        //-- the parameters for the current optimum
887
888        float volBack = volLeft;
889        float volFront = volRight;
890        float newRenderCost = nTotalObjects * totalVol;
891
892#ifdef _DEBUG
893        ofstream sumStats;
894        ofstream vollStats;
895        ofstream volrStats;
896
897        const bool printStats =
898                PrepareOutput(axis, mBvhStats.Leaves(), sumStats, vollStats, volrStats);
899#endif
900
901        ///////////////////////////////////////////////////
902        //-- the sweep heuristics
903        //-- traverse through events and find best split plane
904
905        SortableEntryContainer::const_iterator cit, cit_end = cit_end = mSubdivisionCandidates->end();
906
907        for (cit = mSubdivisionCandidates->begin(); cit != cit_end; ++ cit)
908        {
909                Intersectable *object = (*cit).mObject;
910       
911                // evaluate change in l and r volume
912                // voll = view cells that see only left node (i.e., left pvs)
913                // volr = view cells that see only right node (i.e., right pvs)
914                EvalHeuristicsContribution(object, volLeft, volRight);
915
916                ++ nObjectsLeft;
917                const int nObjectsRight = nTotalObjects - nObjectsLeft;
918
919                // the heuristics
920            const float sum = volLeft * (float)nObjectsLeft +
921                                                  volRight * (float)nObjectsRight;
922
923#ifdef _DEBUG
924                if (printStats)
925                {
926                        PrintHeuristics(nObjectsRight, sum, volLeft, volRight, viewSpaceVol,
927                                                        sumStats, vollStats, volrStats);
928                }
929#endif
930
931                if (sum < newRenderCost)
932                {
933                        newRenderCost = sum;
934
935                        volBack = volLeft;
936                        volFront = volRight;
937
938                        // objects belongs to left side now
939                        for (; backObjectsStart != (cit + 1); ++ backObjectsStart);
940                }
941        }
942
943        ////////////////////////////////////////////
944        //-- assign object to front and back volume
945
946        // belongs to back bv
947        for (cit = mSubdivisionCandidates->begin(); cit != backObjectsStart; ++ cit)
948        {
949                objectsBack.push_back((*cit).mObject);
950        }
951        // belongs to front bv
952        for (cit = backObjectsStart; cit != cit_end; ++ cit)
953        {
954                objectsFront.push_back((*cit).mObject);
955        }
956
957        // render cost of the old parent
958        const float oldRenderCost = (float)nTotalObjects * totalVol + Limits::Small;
959        // the relative cost ratio
960        const float ratio = newRenderCost / oldRenderCost;
961
962#ifdef _DEBUG
963        Debug << "\n§§§§ eval local cost §§§§" << endl
964                  << "back pvs: " << (int)objectsBack.size() << " front pvs: " << (int)objectsFront.size() << " total pvs: " << nTotalObjects << endl
965                  << "back p: " << volBack / viewSpaceVol << " front p " << volFront / viewSpaceVol << " p: " << totalVol / viewSpaceVol << endl
966                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
967                  << "render cost decrease: " << oldRenderCost / viewSpaceVol - newRenderCost / viewSpaceVol << endl;
968#endif
969
970        return ratio;
971}
972
973
974void BvHierarchy::PrepareLocalSubdivisionCandidates(const BvhTraversalData &tData,
975                                                                                                        const int axis)                                                                                 
976{
977        //-- insert object queries
978        ObjectContainer *objects = mUseGlobalSorting ? tData.mSortedObjects[axis] : &tData.mNode->mObjects;
979
980        CreateLocalSubdivisionCandidates(*objects, &mSubdivisionCandidates, !mUseGlobalSorting, axis);
981}
982
983
984void BvHierarchy::CreateLocalSubdivisionCandidates(const ObjectContainer &objects,
985                                                                                                  SortableEntryContainer **subdivisionCandidates,
986                                                                                                  const bool sort,
987                                                                                                  const int axis)
988{
989        (*subdivisionCandidates)->clear();
990
991        // compute requested size and look if subdivision candidate has to be recomputed
992        const int requestedSize = (int)objects.size() * 2;
993       
994        // creates a sorted split candidates array
995        if ((*subdivisionCandidates)->capacity() > 500000 &&
996                requestedSize < (int)((*subdivisionCandidates)->capacity() / 10) )
997        {
998        delete (*subdivisionCandidates);
999                (*subdivisionCandidates) = new SortableEntryContainer;
1000        }
1001
1002        (*subdivisionCandidates)->reserve(requestedSize);
1003
1004        ObjectContainer::const_iterator oit, oit_end = objects.end();
1005
1006        for (oit = objects.begin(); oit < oit_end; ++ oit)
1007        {
1008                Intersectable *object = *oit;
1009                const AxisAlignedBox3 &box = object->GetBox();
1010                const float midPt = (box.Min(axis) + box.Max(axis)) * 0.5f;
1011
1012                (*subdivisionCandidates)->push_back(SortableEntry(object, midPt));
1013        }
1014
1015        if (sort)
1016        {
1017                stable_sort((*subdivisionCandidates)->begin(), (*subdivisionCandidates)->end());
1018        }
1019}
1020
1021
1022const BvhStatistics &BvHierarchy::GetStatistics() const
1023{
1024        return mBvhStats;
1025}
1026
1027
1028float BvHierarchy::PrepareHeuristics(const BvhTraversalData &tData, const int axis)
1029{       
1030        BvhLeaf *leaf = tData.mNode;
1031        float vol = 0;
1032
1033    // sort so we can use a sweep from right to left
1034        PrepareLocalSubdivisionCandidates(tData, axis);
1035       
1036        // collect and mark the view cells as belonging to front pvs
1037        ViewCellContainer viewCells;
1038        CollectViewCells(tData.mNode->mObjects, viewCells, true);
1039                       
1040        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
1041        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
1042        {
1043                vol += (*vit)->GetVolume();
1044        }
1045
1046        // we will mail view cells switching to the back side
1047        ViewCell::NewMail();
1048       
1049        return vol;
1050}
1051///////////////////////////////////////////////////////////
1052
1053
1054void BvHierarchy::EvalHeuristicsContribution(Intersectable *obj,
1055                                                                                         float &volLeft,
1056                                                                                         float &volRight)
1057{
1058        // collect all view cells associated with this objects
1059        // (also multiple times, if they are pierced by several rays)
1060        ViewCellContainer viewCells;
1061        const bool useMailboxing = false;
1062
1063        CollectViewCells(obj, viewCells, useMailboxing);
1064
1065        // classify view cells and compute volume contri accordingly
1066        // possible view cell classifications:
1067        // view cell mailed => view cell can be seen from left child node
1068        // view cell counter > 0 view cell can be seen from right child node
1069        // combined: view cell volume belongs to both nodes
1070        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
1071       
1072        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
1073        {
1074                // view cells can also be seen from left child node
1075                ViewCell *viewCell = *vit;
1076
1077                const float vol = viewCell->GetVolume();
1078
1079                if (!viewCell->Mailed())
1080                {
1081                        viewCell->Mail();
1082                        // we now see view cell from both nodes
1083                        // => add volume to left node
1084                        volLeft += vol;
1085                }
1086
1087                // last reference into the right node
1088                if (-- viewCell->mCounter == 0)
1089                {       
1090                        // view cell was previously seen from both nodes  =>
1091                        // remove volume from right node
1092                        volRight -= vol;
1093                }
1094        }
1095}
1096
1097
1098void BvHierarchy::SetViewCellsManager(ViewCellsManager *vcm)
1099{
1100        mViewCellsManager = vcm;
1101}
1102
1103
1104AxisAlignedBox3 BvHierarchy::GetBoundingBox() const
1105{
1106        return mBoundingBox;
1107}
1108
1109
1110float BvHierarchy::SelectObjectPartition(const BvhTraversalData &tData,
1111                                                                                 ObjectContainer &frontObjects,
1112                                                                                 ObjectContainer &backObjects)
1113{
1114        ObjectContainer nFrontObjects[3];
1115        ObjectContainer nBackObjects[3];
1116        float nCostRatio[3];
1117
1118        int sAxis = 0;
1119        int bestAxis = -1;
1120
1121        if (mOnlyDrivingAxis)
1122        {
1123                const AxisAlignedBox3 box = tData.mNode->GetBoundingBox();
1124                sAxis = box.Size().DrivingAxis();
1125        }
1126       
1127        ////////////////////////////////////////////
1128        //-- evaluate split cost for all three axis
1129       
1130        for (int axis = 0; axis < 3; ++ axis)
1131        {
1132                if (!mOnlyDrivingAxis || (axis == sAxis))
1133                {
1134                        if (mUseCostHeuristics)
1135                        {
1136                                //////////////////////////////////
1137                //-- split objects using heuristics
1138                               
1139                                if (mHierarchyManager->GetViewSpaceSubdivisionType() ==
1140                                        HierarchyManager::KD_BASED_VIEWSPACE_SUBDIV)
1141                                {
1142                                        //-- heuristics using objects weighted by view cells volume
1143                                        nCostRatio[axis] =
1144                                                EvalLocalCostHeuristics(
1145                                                                                                tData,
1146                                                                                                axis,
1147                                                                                                nFrontObjects[axis],
1148                                                                                                nBackObjects[axis]);
1149                                }
1150                                else
1151                                {
1152                                        //-- use surface area heuristic because view cells not constructed yet                         
1153                                        nCostRatio[axis] =
1154                                                EvalSah(
1155                                                                tData,
1156                                                                axis,
1157                                                                nFrontObjects[axis],
1158                                                                nBackObjects[axis]);
1159                                }
1160                        }
1161                        else
1162                        {
1163                                //-- split objects using some simple criteria
1164                                nCostRatio[axis] =
1165                                        EvalLocalObjectPartition(
1166                                                                                         tData,
1167                                                                                         axis,
1168                                                                                         nFrontObjects[axis],
1169                                                                                         nBackObjects[axis]);
1170                        }
1171
1172                        if (bestAxis == -1)
1173                        {
1174                                bestAxis = axis;
1175                        }
1176                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1177                        {
1178                                bestAxis = axis;
1179                        }
1180                }
1181        }
1182
1183        /////////////////////////////////////
1184        //-- assign values
1185
1186        frontObjects = nFrontObjects[bestAxis];
1187        backObjects = nBackObjects[bestAxis];
1188
1189        //Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
1190        return nCostRatio[bestAxis];
1191}
1192
1193
1194int BvHierarchy::AssociateObjectsWithRays(const VssRayContainer &rays) const
1195{
1196        int nRays = 0;
1197        VssRayContainer::const_iterator rit, rit_end = rays.end();
1198
1199        VssRay::NewMail();
1200
1201    for (rit = rays.begin(); rit != rays.end(); ++ rit)
1202        {
1203                VssRay *ray = (*rit);
1204
1205                if (ray->mTerminationObject)
1206                {
1207                        ray->mTerminationObject->mVssRays.push_back(ray);
1208                        if (!ray->Mailed())
1209                        {
1210                                ray->Mail();
1211                                ++ nRays;
1212                        }
1213                }
1214
1215                if (1 && ray->mOriginObject)
1216                {
1217                        ray->mOriginObject->mVssRays.push_back(ray);
1218
1219                        if (!ray->Mailed())
1220                        {
1221                                ray->Mail();
1222                                ++ nRays;
1223                        }
1224                }
1225        }
1226
1227        return nRays;
1228}
1229
1230
1231void BvHierarchy::PrintSubdivisionStats(const SubdivisionCandidate &sc)
1232{
1233        const float costDecr = sc.GetRenderCostDecrease();     
1234
1235        mSubdivisionStats
1236                        << "#Leaves\n" << mBvhStats.Leaves()
1237                        << "#RenderCostDecrease\n" << costDecr << endl
1238                        << "#TotalRenderCost\n" << mTotalCost << endl;
1239                        //<< "#AvgRenderCost\n" << avgRenderCost << endl;
1240}
1241
1242
1243void BvHierarchy::CollectRays(const ObjectContainer &objects,
1244                                                          VssRayContainer &rays) const
1245{
1246        VssRay::NewMail();
1247        ObjectContainer::const_iterator oit, oit_end = objects.end();
1248
1249        // evaluate reverse pvs and view cell volume on left and right cell
1250        // note: should I take all leaf objects or rather the objects hit by rays?
1251        for (oit = objects.begin(); oit != oit_end; ++ oit)
1252        {
1253                Intersectable *obj = *oit;
1254                VssRayContainer::const_iterator rit, rit_end = obj->mVssRays.end();
1255
1256                for (rit = obj->mVssRays.begin(); rit < rit_end; ++ rit)
1257                {
1258                        VssRay *ray = (*rit);
1259
1260                        if (!ray->Mailed())
1261                        {
1262                                ray->Mail();
1263                                rays.push_back(ray);
1264                        }
1265                }
1266        }
1267}
1268
1269
1270float BvHierarchy::EvalRenderCost(const ObjectContainer &objects) const
1271{
1272        if (mHierarchyManager->GetViewSpaceSubdivisionType() == HierarchyManager::NO_VIEWSPACE_SUBDIV)
1273        {
1274                if (objects.empty())
1275                        return 0.0f;
1276
1277                /////////////////////////////
1278                //-- surface area heuristics
1279
1280                const AxisAlignedBox3 box = EvalBoundingBox(objects);
1281                const float area = box.SurfaceArea();
1282
1283                return (float)objects.size() * area;
1284        }
1285        else
1286        {       /////////////////////////////
1287                //-- render cost heuristics
1288
1289                const float viewSpaceVol = mHierarchyManager->GetViewSpaceBox().GetVolume();
1290                // probability that view point lies in a view cell which sees this node
1291                const float p = EvalViewCellsVolume(objects) / viewSpaceVol;   
1292
1293                return (float)objects.size() * p;
1294        }
1295}
1296
1297
1298AxisAlignedBox3 BvHierarchy::EvalBoundingBox(const ObjectContainer &objects,
1299                                                                                         const AxisAlignedBox3 *parentBox) const
1300{
1301        // if there are no objects in this box, box size is set to parent box size.
1302        // Question: Invalidate box instead?
1303        if (parentBox && objects.empty())
1304                return *parentBox;
1305
1306        AxisAlignedBox3 box;
1307        box.Initialize();
1308
1309        ObjectContainer::const_iterator oit, oit_end = objects.end();
1310
1311        for (oit = objects.begin(); oit != oit_end; ++ oit)
1312        {
1313                Intersectable *obj = *oit;
1314                // grow bounding box to include all objects
1315                box.Include(obj->GetBox());
1316        }
1317
1318        return box;
1319}
1320
1321
1322void BvHierarchy::CollectLeaves(vector<BvhLeaf *> &leaves) const
1323{
1324        stack<BvhNode *> nodeStack;
1325        nodeStack.push(mRoot);
1326
1327        while (!nodeStack.empty())
1328        {
1329                BvhNode *node = nodeStack.top();
1330                nodeStack.pop();
1331
1332                if (node->IsLeaf())
1333                {
1334                        BvhLeaf *leaf = (BvhLeaf *)node;
1335                        leaves.push_back(leaf);
1336                }
1337                else
1338                {
1339                        BvhInterior *interior = (BvhInterior *)node;
1340
1341                        nodeStack.push(interior->GetBack());
1342                        nodeStack.push(interior->GetFront());
1343                }
1344        }
1345}
1346
1347
1348AxisAlignedBox3 BvHierarchy::GetBoundingBox(BvhNode *node) const
1349{
1350        return node->GetBoundingBox();
1351}
1352
1353
1354void BvHierarchy::CollectViewCells(const ObjectContainer &objects,
1355                                                                   ViewCellContainer &viewCells,
1356                                                                   const bool setCounter) const
1357{
1358        // no view cells yet
1359        if (mHierarchyManager->GetViewSpaceSubdivisionType() ==
1360                HierarchyManager::NO_VIEWSPACE_SUBDIV)
1361                return;
1362
1363        ViewCell::NewMail();
1364        ObjectContainer::const_iterator oit, oit_end = objects.end();
1365
1366        // loop through all object and collect view cell pvs of this node
1367        for (oit = objects.begin(); oit != oit_end; ++ oit)
1368        {
1369                CollectViewCells(*oit, viewCells, true, setCounter);
1370        }
1371}
1372
1373
1374void BvHierarchy::CollectViewCells(Intersectable *obj,
1375                                                                   ViewCellContainer &viewCells,
1376                                                                   const bool useMailBoxing,
1377                                                                   const bool setCounter) const
1378{
1379        VssRayContainer::const_iterator rit, rit_end = obj->mVssRays.end();
1380
1381        for (rit = obj->mVssRays.begin(); rit < rit_end; ++ rit)
1382        {
1383                VssRay *ray = (*rit);
1384                ViewCellContainer tmpViewCells;
1385       
1386                mHierarchyManager->mVspTree->GetViewCells(*ray, tmpViewCells);
1387
1388                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
1389
1390                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
1391                {
1392                        VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);
1393
1394                        // store view cells
1395                        if (!useMailBoxing || !vc->Mailed())
1396                        {
1397                                if (useMailBoxing)
1398                                {
1399                                        vc->Mail();
1400                                        if (setCounter)
1401                                        {
1402                                                vc->mCounter = 0;
1403                                        }
1404                                }
1405                                viewCells.push_back(vc);
1406                        }
1407                       
1408                        if (setCounter)
1409                        {
1410                                ++ vc->mCounter;
1411                        }
1412                }
1413        }
1414}
1415
1416
1417void BvHierarchy::CollectDirtyCandidates(BvhSubdivisionCandidate *sc,
1418                                                                                 vector<SubdivisionCandidate *> &dirtyList)
1419{
1420        BvhTraversalData &tData = sc->mParentData;
1421        BvhLeaf *node = tData.mNode;
1422       
1423        ViewCellContainer viewCells;
1424        CollectViewCells(node->mObjects, viewCells);
1425
1426        // split candidates handling
1427        // these view cells  are thrown into dirty list
1428        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
1429
1430        Debug << "collecting " << (int)viewCells.size() << " dirty candidates" << endl;
1431
1432        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
1433        {
1434                VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);
1435                VspLeaf *leaf = vc->mLeaf;
1436                SubdivisionCandidate *candidate = leaf->GetSubdivisionCandidate();
1437               
1438                if (candidate) // is this leaf still a split candidate?
1439                {
1440                        dirtyList.push_back(candidate);
1441                }
1442        }
1443}
1444
1445
1446BvhNode *BvHierarchy::GetRoot() const
1447{
1448        return mRoot;
1449}
1450
1451
1452bool BvHierarchy::IsObjectInLeaf(BvhLeaf *leaf, Intersectable *object) const
1453{
1454        ObjectContainer::const_iterator oit =
1455                lower_bound(leaf->mObjects.begin(), leaf->mObjects.end(), object, ilt);
1456                               
1457        // objects sorted by id
1458        if ((oit != leaf->mObjects.end()) && ((*oit)->GetId() == object->GetId()))
1459        {
1460                return true;
1461        }
1462        else
1463        {
1464                return false;
1465        }
1466}
1467
1468
1469BvhLeaf *BvHierarchy::GetLeaf(Intersectable *object, BvhNode *node) const
1470{
1471        // rather use the simple version
1472        if (!object) return NULL;
1473        return object->mBvhLeaf;
1474       
1475        ///////////////////////////////////////
1476        // start from root of tree
1477       
1478        if (node == NULL)
1479                node = mRoot;
1480       
1481        vector<BvhLeaf *> leaves;
1482
1483        stack<BvhNode *> nodeStack;
1484        nodeStack.push(node);
1485 
1486        BvhLeaf *leaf = NULL;
1487 
1488        while (!nodeStack.empty()) 
1489        {
1490                BvhNode *node = nodeStack.top();
1491                nodeStack.pop();
1492       
1493                if (node->IsLeaf())
1494                {
1495                        leaf = dynamic_cast<BvhLeaf *>(node);
1496
1497                        if (IsObjectInLeaf(leaf, object))
1498                        {
1499                                return leaf;
1500                        }
1501                }
1502                else   
1503                {       
1504                        // find point
1505                        BvhInterior *interior = dynamic_cast<BvhInterior *>(node);
1506       
1507                        if (interior->GetBack()->GetBoundingBox().Includes(object->GetBox()))
1508                        {
1509                                nodeStack.push(interior->GetBack());
1510                        }
1511                       
1512                        // search both sides as we are using bounding volumes
1513                        if (interior->GetFront()->GetBoundingBox().Includes(object->GetBox()))
1514                        {
1515                                nodeStack.push(interior->GetFront());
1516                        }
1517                }
1518        }
1519 
1520        return leaf;
1521}
1522
1523
1524BvhIntersectable *BvHierarchy::GetOrCreateBvhIntersectable(BvhNode *node)
1525{
1526        // search nodes
1527        std::map<BvhNode *, BvhIntersectable *>::
1528                const_iterator it = mBvhIntersectables.find(node);
1529
1530        if (it != mBvhIntersectables.end())
1531        {
1532                return (*it).second;
1533        }
1534
1535        // not in map => create new entry
1536        BvhIntersectable *bvhObj = new BvhIntersectable(node);
1537        mBvhIntersectables[node] = bvhObj;
1538
1539        return bvhObj;
1540}
1541
1542
1543bool BvHierarchy::Export(OUT_STREAM &stream)
1544{
1545        ExportNode(mRoot, stream);
1546
1547        return true;
1548}
1549
1550
1551void BvHierarchy::ExportObjects(BvhLeaf *leaf, OUT_STREAM &stream)
1552{
1553        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
1554        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
1555        {
1556                stream << (*oit)->GetId() << " ";
1557        }
1558}
1559
1560
1561void BvHierarchy::ExportNode(BvhNode *node, OUT_STREAM &stream)
1562{
1563        if (node->IsLeaf())
1564        {
1565                BvhLeaf *leaf = dynamic_cast<BvhLeaf *>(node);
1566                const AxisAlignedBox3 box = leaf->GetBoundingBox();
1567                stream << "<Leaf"
1568                           << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
1569                           << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\""
1570                           << " objects=\"";
1571               
1572                //-- export objects
1573                ExportObjects(leaf, stream);
1574               
1575                stream << "\" />" << endl;
1576        }
1577        else
1578        {       
1579                BvhInterior *interior = dynamic_cast<BvhInterior *>(node);
1580                const AxisAlignedBox3 box = interior->GetBoundingBox();
1581
1582                stream << "<Interior"
1583                           << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
1584                           << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z
1585                           << "\">" << endl;
1586
1587                ExportNode(interior->GetBack(), stream);
1588                ExportNode(interior->GetFront(), stream);
1589
1590                stream << "</Interior>" << endl;
1591        }
1592}
1593
1594
1595float BvHierarchy::EvalViewCellsVolume(const ObjectContainer &objects) const
1596{
1597        float vol = 0;
1598
1599        ViewCellContainer viewCells;
1600        CollectViewCells(objects, viewCells);
1601
1602        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
1603
1604        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
1605        {
1606                vol += (*vit)->GetVolume();
1607        }
1608
1609        return vol;
1610}
1611
1612
1613void BvHierarchy::CreateRoot(const ObjectContainer &objects)
1614{
1615        //-- create new root
1616        AxisAlignedBox3 box = EvalBoundingBox(objects);
1617        BvhLeaf *bvhleaf = new BvhLeaf(box, NULL, (int)objects.size());
1618        bvhleaf->mObjects = objects;
1619        mRoot = bvhleaf;
1620
1621        // associate root with current objects
1622        AssociateObjectsWithLeaf(bvhleaf);
1623}
1624
1625/*
1626Mesh *BvHierarchy::MergeLeafToMesh()
1627void BvHierarchy::MergeLeavesToMeshes()
1628{
1629        vector<BvhLeaf *> leaves;
1630        CollectLeaves(leaves);
1631
1632        vector<BvhLeaf *>::const_iterator lit, lit_end = leaves.end();
1633
1634        for (lit = leaves.begin(); lit != lit_end; ++ lit)
1635        {
1636                Mesh *mesh = MergeLeafToMesh(*lit);
1637        }
1638}*/
1639
1640
1641SubdivisionCandidate *BvHierarchy::PrepareConstruction(const VssRayContainer &sampleRays,
1642                                                                                                           const ObjectContainer &objects)
1643{
1644        ///////////////////////////////////////////////////////////////
1645        //-- note matt: we assume that we have objects sorted by their id =>
1646        //-- we don't have to sort them here and an binary search
1647        //-- for identifying if a object is in a leaf.
1648
1649        mBvhStats.Reset();
1650        mBvhStats.Start();
1651        mBvhStats.nodes = 1;
1652
1653        // store pointer to this tree
1654        BvhSubdivisionCandidate::sBvHierarchy = this;
1655        mBvhStats.nodes = 1;
1656        mGlobalCostMisses = 0;
1657
1658        // compute bounding box from objects
1659        // note: we assume that root was already created
1660        mBoundingBox = mRoot->GetBoundingBox();
1661        BvhLeaf *bvhleaf = dynamic_cast<BvhLeaf *>(mRoot);
1662
1663        // multiply termination criterium for comparison,
1664        // so it can be set between zero and one and
1665        // no division is necessary during traversal
1666
1667#if PROBABILIY_IS_BV_VOLUME
1668        mTermMinProbability *= mBoundingBox.GetVolume();
1669        // probability that bounding volume is seen
1670        const float prop = GetBoundingBox().GetVolume();
1671#else
1672        mTermMinProbability *= mVspTree->GetBoundingBox().GetVolume();
1673        // probability that volume is "seen" from the view cells
1674        const float prop = EvalViewCellsVolume(objects);
1675#endif
1676
1677        // only rays intersecting objects in node are interesting
1678        const int nRays = AssociateObjectsWithRays(sampleRays);
1679        //Debug << "using " << nRays << " of " << (int)sampleRays.size() << " rays" << endl;
1680
1681        // create bvh traversal data
1682        BvhTraversalData oData(bvhleaf, 0, prop, nRays);
1683
1684        // create sorted object lists for the first data
1685        if (mUseGlobalSorting)
1686        {
1687                CreateInitialSortedObjectList(oData);
1688        }
1689       
1690
1691        ////////////////////////////////////////////////////
1692        //-- add first candidate for object space partition     
1693
1694        BvhSubdivisionCandidate *oSubdivisionCandidate =
1695                new BvhSubdivisionCandidate(oData);
1696
1697        EvalSubdivisionCandidate(*oSubdivisionCandidate);
1698        bvhleaf->SetSubdivisionCandidate(oSubdivisionCandidate);
1699
1700        const float viewSpaceVol = mHierarchyManager->GetViewSpaceBox().GetVolume();
1701        mTotalCost = (float)objects.size() * prop / viewSpaceVol;
1702
1703        PrintSubdivisionStats(*oSubdivisionCandidate);
1704
1705        return oSubdivisionCandidate;
1706}
1707
1708
1709void BvHierarchy::CreateInitialSortedObjectList(BvhTraversalData &tData)
1710{
1711        SortableEntryContainer *sortedObjects = new SortableEntryContainer();
1712
1713        // we sort the objects as a preprocess so they don't have
1714        // to be sorted for each split
1715        for (int i = 0; i < 3; ++ i)
1716        {
1717                CreateLocalSubdivisionCandidates(tData.mNode->mObjects, &sortedObjects, true, i);
1718                       
1719                tData.mSortedObjects[i] = new ObjectContainer();
1720                tData.mSortedObjects[i]->reserve((int)sortedObjects->size());
1721
1722                SortableEntryContainer::const_iterator oit, oit_end = sortedObjects->end();
1723                for (oit = sortedObjects->begin(); oit != oit_end; ++ oit)
1724                {
1725                        tData.mSortedObjects[i]->push_back((*oit).mObject);
1726                }
1727        }
1728
1729        delete sortedObjects;
1730}
1731
1732
1733void BvHierarchy::AssignSortedObjects(const BvhSubdivisionCandidate &sc,
1734                                                                          BvhTraversalData &frontData,
1735                                                                          BvhTraversalData &backData)
1736{
1737        Intersectable::NewMail();
1738
1739        // we sorted the objects as a preprocess so they don't have
1740        // to be sorted for each split
1741        ObjectContainer::const_iterator fit, fit_end = sc.mFrontObjects.end();
1742
1743        for (fit = sc.mFrontObjects.begin(); fit != fit_end; ++ fit)
1744        {
1745                (*fit)->Mail();
1746        }
1747
1748        for (int i = 0; i < 3; ++ i)
1749        {
1750                frontData.mSortedObjects[i] = new ObjectContainer();
1751                backData.mSortedObjects[i] = new ObjectContainer();
1752
1753                frontData.mSortedObjects[i]->reserve((int)sc.mFrontObjects.size());
1754                backData.mSortedObjects[i]->reserve((int)sc.mFrontObjects.size());
1755
1756                ObjectContainer::const_iterator oit, oit_end = sc.mParentData.mSortedObjects[i]->end();
1757
1758                for (oit = sc.mParentData.mSortedObjects[i]->begin(); oit != oit_end; ++ oit)
1759                {
1760                        if ((*oit)->Mailed())
1761                        {
1762                                frontData.mSortedObjects[i]->push_back(*oit);
1763                        }
1764                        else
1765                        {
1766                                backData.mSortedObjects[i]->push_back(*oit);
1767                        }
1768                }
1769        }
1770}
1771
1772
1773void BvhStatistics::Print(ostream &app) const
1774{
1775        app << "=========== BvHierarchy statistics ===============\n";
1776
1777        app << setprecision(4);
1778
1779        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
1780
1781        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
1782
1783        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
1784
1785        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
1786
1787        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits << endl;
1788
1789        app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
1790                << maxCostNodes * 100 / (double)Leaves() << endl;
1791
1792        app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
1793                << minProbabilityNodes * 100 / (double)Leaves() << endl;
1794
1795
1796        //////////////////////////////////////////////////
1797       
1798        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
1799                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
1800       
1801        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
1802
1803        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
1804
1805        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
1806
1807       
1808        ////////////////////////////////////////////////////////
1809       
1810        app << "#N_PMINOBJECTSLEAVES  ( Percentage of leaves with mininum objects )\n"
1811                << minObjectsNodes * 100 / (double)Leaves() << endl;
1812
1813        app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" << maxObjectRefs << "\n";
1814
1815        app << "#N_MINOBJECTREFS  ( Min number of object refs / leaf )\n" << minObjectRefs << "\n";
1816
1817        app << "#N_EMPTYLEAFS ( Empty leafs )\n" << emptyNodes << "\n";
1818       
1819        app << "#N_PAVGOBJECTSLEAVES  ( average object refs / leaf)\n" << AvgObjectRefs() << endl;
1820
1821
1822        ////////////////////////////////////////////////////////
1823       
1824        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with mininum rays )\n"
1825                << minRaysNodes * 100 / (double)Leaves() << endl;
1826
1827        app << "#N_MAXRAYREFS  ( Max number of ray refs / leaf )\n" << maxRayRefs << "\n";
1828
1829        app << "#N_MINRAYREFS  ( Min number of ray refs / leaf )\n" << minRayRefs << "\n";
1830       
1831        app << "#N_PAVGRAYLEAVES  ( average ray refs / leaf )\n" << AvgRayRefs() << endl;
1832       
1833        app << "#N_PAVGRAYCONTRIBLEAVES  ( Average ray contribution)\n" <<
1834                rayRefs / (double)objectRefs << endl;
1835
1836        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
1837                maxRayContriNodes * 100 / (double)Leaves() << endl;
1838
1839        app << "========== END OF BvHierarchy statistics ==========\n";
1840}
1841
1842
1843}
Note: See TracBrowser for help on using the repository browser.