source: GTP/trunk/Lib/Vis/Preprocessing/src/VspOspTree.cpp @ 1012

Revision 1012, 51.3 KB checked in by mattausch, 18 years ago (diff)
Line 
1#include <stack>
2#include <time.h>
3#include <iomanip>
4
5#include "ViewCell.h"
6#include "Plane3.h"
7#include "VspOspTree.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
19namespace GtpVisibilityPreprocessor {
20
21#define USE_FIXEDPOINT_T 0
22
23
24//-- static members
25
26int VspOspTree::sFrontId = 0;
27int VspOspTree::sBackId = 0;
28int VspOspTree::sFrontAndBackId = 0;
29
30
31
32// pvs penalty can be different from pvs size
33inline static float EvalPvsPenalty(const int pvs,
34                                                                   const int lower,
35                                                                   const int upper)
36{
37        // clamp to minmax values
38        if (pvs < lower)
39                return (float)lower;
40        if (pvs > upper)
41                return (float)upper;
42
43        return (float)pvs;
44}
45
46
47int VspNode::sMailId = 1;
48
49
50/******************************************************************/
51/*                  class VspNode implementation                  */
52/******************************************************************/
53
54
55VspNode::VspNode():
56mParent(NULL), mTreeValid(true), mTimeStamp(0)
57{}
58
59
60VspNode::VspNode(VspInterior *parent):
61mParent(parent), mTreeValid(true)
62{}
63
64
65bool VspNode::IsRoot() const
66{
67        return mParent == NULL;
68}
69
70
71VspInterior *VspNode::GetParent()
72{
73        return mParent;
74}
75
76
77void VspNode::SetParent(VspInterior *parent)
78{
79        mParent = parent;
80}
81
82
83bool VspNode::IsSibling(VspNode *n) const
84{
85        return  ((this != n) && mParent &&
86                         (mParent->GetFront() == n) || (mParent->GetBack() == n));
87}
88
89
90int VspNode::GetDepth() const
91{
92        int depth = 0;
93        VspNode *p = mParent;
94       
95        while (p)
96        {
97                p = p->mParent;
98                ++ depth;
99        }
100
101        return depth;
102}
103
104
105bool VspNode::TreeValid() const
106{
107        return mTreeValid;
108}
109
110
111void VspNode::SetTreeValid(const bool v)
112{
113        mTreeValid = v;
114}
115
116
117/****************************************************************/
118/*              class VspInterior implementation                */
119/****************************************************************/
120
121
122VspInterior::VspInterior(const AxisAlignedPlane &plane):
123mPlane(plane), mFront(NULL), mBack(NULL)
124{}
125
126
127VspInterior::~VspInterior()
128{
129        DEL_PTR(mFront);
130        DEL_PTR(mBack);
131}
132
133
134bool VspInterior::IsLeaf() const
135{
136        return false;
137}
138
139
140VspNode *VspInterior::GetBack()
141{
142        return mBack;
143}
144
145
146VspNode *VspInterior::GetFront()
147{
148        return mFront;
149}
150
151
152AxisAlignedPlane VspInterior::GetPlane() const
153{
154        return mPlane;
155}
156
157
158float VspInterior::GetPosition() const
159{
160        return mPlane.mPosition;
161}
162
163
164int VspInterior::GetAxis() const
165{
166        return mPlane.mAxis;
167}
168
169
170void VspInterior::ReplaceChildLink(VspNode *oldChild, VspNode *newChild)
171{
172        if (mBack == oldChild)
173                mBack = newChild;
174        else
175                mFront = newChild;
176}
177
178
179void VspInterior::SetupChildLinks(VspNode *b, VspNode *f)
180{
181    mBack = b;
182    mFront = f;
183}
184
185
186AxisAlignedBox3 VspInterior::GetBox() const
187{
188        return mBox;
189}
190
191
192void VspInterior::SetBox(const AxisAlignedBox3 &box)
193{
194        mBox = box;
195}
196
197
198int VspInterior::Type() const
199{
200        return Interior;
201}
202
203/****************************************************************/
204/*                  class VspLeaf implementation                */
205/****************************************************************/
206
207
208VspLeaf::VspLeaf(): mViewCell(NULL), mPvs(NULL)
209{
210}
211
212
213VspLeaf::~VspLeaf()
214{
215        DEL_PTR(mPvs);
216}
217
218int VspLeaf::Type() const
219{
220        return Leaf;
221}
222
223
224
225VspLeaf::VspLeaf(ViewCellLeaf *viewCell):
226mViewCell(viewCell)
227{
228}
229
230
231VspLeaf::VspLeaf(VspInterior *parent):
232VspNode(parent), mViewCell(NULL), mPvs(NULL)
233{}
234
235
236
237VspLeaf::VspLeaf(VspInterior *parent, ViewCellLeaf *viewCell):
238VspNode(parent), mViewCell(viewCell), mPvs(NULL)
239{
240}
241
242ViewCellLeaf *VspLeaf::GetViewCell() const
243{
244        return mViewCell;
245}
246
247void VspLeaf::SetViewCell(ViewCellLeaf *viewCell)
248{
249        mViewCell = viewCell;
250}
251
252
253bool VspLeaf::IsLeaf() const
254{
255        return true;
256}
257
258
259/******************************************************************************/
260/*                       class VspOspTree implementation                      */
261/******************************************************************************/
262
263
264VspOspTree::VspOspTree():
265mRoot(NULL),
266mOutOfBoundsCell(NULL),
267mStoreRays(false),
268mUseRandomAxis(false),
269mTimeStamp(1)
270{
271        bool randomize = false;
272        Environment::GetSingleton()->GetBoolValue("VspOspTree.Construction.randomize", randomize);
273        if (randomize)
274                Randomize(); // initialise random generator for heuristics
275
276        //-- termination criteria for autopartition
277        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.maxDepth", mTermMaxDepth);
278        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.minPvs", mTermMinPvs);
279        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.minRays", mTermMinRays);
280        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.minProbability", mTermMinProbability);
281        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxRayContribution", mTermMaxRayContribution);
282        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxCostRatio", mTermMaxCostRatio);
283        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.missTolerance", mTermMissTolerance);
284        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.maxViewCells", mMaxViewCells);
285
286        //-- max cost ratio for early tree termination
287        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.maxCostRatio", mTermMaxCostRatio);
288
289        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
290        Environment::GetSingleton()->GetIntValue("VspOspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
291
292        // HACK//mTermMinPolygons = 25;
293
294        //-- factors for bsp tree split plane heuristics
295        Environment::GetSingleton()->GetFloatValue("VspOspTree.Termination.ct_div_ci", mCtDivCi);
296
297        //-- partition criteria
298        Environment::GetSingleton()->GetFloatValue("VspOspTree.Construction.epsilon", mEpsilon);
299
300        // if only the driving axis is used for axis aligned split
301        Environment::GetSingleton()->GetBoolValue("VspOspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
302       
303        //Environment::GetSingleton()->GetFloatValue("VspOspTree.maxTotalMemory", mMaxTotalMemory);
304        Environment::GetSingleton()->GetFloatValue("VspOspTree.maxStaticMemory", mMaxMemory);
305
306        Environment::GetSingleton()->GetBoolValue("VspOspTree.useCostHeuristics", mUseCostHeuristics);
307        Environment::GetSingleton()->GetBoolValue("VspOspTree.simulateOctree", mCirculatingAxis);
308        Environment::GetSingleton()->GetBoolValue("VspOspTree.useRandomAxis", mUseRandomAxis);
309        Environment::GetSingleton()->GetIntValue("VspOspTree.nodePriorityQueueType", mNodePriorityQueueType);
310       
311        char subdivisionStatsLog[100];
312        Environment::GetSingleton()->GetStringValue("VspOspTree.subdivisionStats", subdivisionStatsLog);
313        mSubdivisionStats.open(subdivisionStatsLog);
314
315        Environment::GetSingleton()->GetFloatValue("VspOspTree.Construction.minBand", mMinBand);
316        Environment::GetSingleton()->GetFloatValue("VspOspTree.Construction.maxBand", mMaxBand);
317       
318
319        //-- debug output
320
321        Debug << "******* VSP BSP options ******** " << endl;
322    Debug << "max depth: " << mTermMaxDepth << endl;
323        Debug << "min PVS: " << mTermMinPvs << endl;
324        Debug << "min probabiliy: " << mTermMinProbability << endl;
325        Debug << "min rays: " << mTermMinRays << endl;
326        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
327        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
328        Debug << "miss tolerance: " << mTermMissTolerance << endl;
329        Debug << "max view cells: " << mMaxViewCells << endl;
330        Debug << "randomize: " << randomize << endl;
331
332        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
333        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
334        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
335        Debug << "max memory: " << mMaxMemory << endl;
336        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
337        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
338        Debug << "use random axis: " << mUseRandomAxis << endl;
339        Debug << "priority queue type: " << mNodePriorityQueueType << endl;
340        Debug << "circulating axis: " << mCirculatingAxis << endl;
341        Debug << "minband: " << mMinBand << endl;
342        Debug << "maxband: " << mMaxBand << endl;
343       
344
345        mSplitCandidates = new vector<SortableEntry>;
346
347        Debug << endl;
348}
349
350
351VspViewCell *VspOspTree::GetOutOfBoundsCell()
352{
353        return mOutOfBoundsCell;
354}
355
356
357VspViewCell *VspOspTree::GetOrCreateOutOfBoundsCell()
358{
359        if (!mOutOfBoundsCell)
360        {
361                mOutOfBoundsCell = new VspViewCell();
362                mOutOfBoundsCell->SetId(-1);
363                mOutOfBoundsCell->SetValid(false);
364        }
365
366        return mOutOfBoundsCell;
367}
368
369
370const VspTreeStatistics &VspOspTree::GetStatistics() const
371{
372        return mVspStats;
373}
374
375
376VspOspTree::~VspOspTree()
377{
378        DEL_PTR(mRoot);
379        DEL_PTR(mSplitCandidates);
380}
381
382
383void VspOspTree::Construct(const VssRayContainer &sampleRays,
384                                                   AxisAlignedBox3 *forcedBoundingBox)
385{
386        mVspStats.nodes = 1;
387        mBox.Initialize();      // initialise BSP tree bounding box
388
389        if (forcedBoundingBox)
390                mBox = *forcedBoundingBox;
391       
392        PolygonContainer polys;
393        RayInfoContainer *rays = new RayInfoContainer();
394
395        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
396
397        long startTime = GetTime();
398
399        cout << "Extracting polygons from rays ... ";
400
401        Intersectable::NewMail();
402
403        int numObj = 0;
404
405        //-- store rays
406        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
407        {
408                VssRay *ray = *rit;
409
410                float minT, maxT;
411
412                static Ray hray;
413                hray.Init(*ray);
414
415                // TODO: not very efficient to implictly cast between rays types
416                if (mBox.GetRaySegment(hray, minT, maxT))
417                {
418                        float len = ray->Length();
419
420                        if (!len)
421                                len = Limits::Small;
422
423                        rays->push_back(RayInfo(ray, minT / len, maxT / len));
424                }
425        }
426
427        mTermMinProbability *= mBox.GetVolume();
428
429        mGlobalCostMisses = 0;
430
431        cout << "finished" << endl;
432
433        // use split cost priority queue
434        Construct(rays);
435}
436
437
438// TODO: return memory usage in MB
439float VspOspTree::GetMemUsage() const
440{
441        return (float)
442                 (sizeof(VspOspTree) +
443                  mVspStats.Leaves() * sizeof(VspLeaf) +
444                  mCreatedViewCells * sizeof(VspViewCell) +
445                  mVspStats.pvs * sizeof(ObjectPvsData) +
446                  mVspStats.Interior() * sizeof(VspInterior) +
447                  mVspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
448}
449
450
451void VspOspTree::Construct(RayInfoContainer *rays)
452{
453        VspOspSplitQueue tQueue;
454
455        mRoot = new VspLeaf();
456
457        const float prop = mBox.GetVolume();
458
459        VspOspTraversalData tData(mRoot,
460                                                          0,
461                                                          rays,
462                              ComputePvsSize(*rays),
463                                                          prop,
464                                                          mBox);
465
466
467        // compute first split candidate
468        VspOspSplitCandidate splitCandidate;
469        EvalSplitCandidate(tData, splitCandidate);
470
471        tQueue.push(splitCandidate);
472
473        mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
474        mTotalPvsSize = tData.mPvs;
475       
476        mSubdivisionStats
477                        << "#ViewCells\n1\n" <<  endl
478                        << "#RenderCostDecrease\n0\n" << endl
479                        << "#SplitCandidateCost\n0\n" << endl
480                        << "#TotalRenderCost\n" << mTotalCost << endl
481                        << "#AvgRenderCost\n" << mTotalPvsSize << endl;
482
483        Debug << "total cost: " << mTotalCost << endl;
484       
485       
486        mVspStats.Start();
487        cout << "Constructing vsp bsp tree ... \n";
488
489        long startTime = GetTime();     
490        int nLeaves = 500;
491        int nViewCells = 500;
492
493        // used for intermediate time measurements and progress
494        long interTime = GetTime();     
495
496        mOutOfMemory = false;
497
498        mCreatedViewCells = 0;
499       
500        while (!tQueue.empty())
501        {
502                splitCandidate = tQueue.top();
503            tQueue.pop();               
504
505                // cost ratio of cost decrease / totalCost
506                float costRatio = splitCandidate.GetCost() / mTotalCost;
507
508                //Debug << "cost ratio: " << costRatio << endl;
509
510                if (costRatio < mTermMinGlobalCostRatio)
511                        ++ mGlobalCostMisses;
512               
513                if (0 && !mOutOfMemory)
514                {
515                        float mem = GetMemUsage();
516
517                        if (mem > mMaxMemory)
518                        {
519                                mOutOfMemory = true;
520                                Debug << "memory limit reached: " << mem << endl;
521                        }
522                }
523
524                // subdivide leaf node
525                VspNode *r = Subdivide(tQueue, splitCandidate);
526
527                if (r == mRoot)
528                        Debug << "VSP BSP tree construction time spent at root: "
529                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
530
531                if (mVspStats.Leaves() >= nLeaves)
532                {
533                        nLeaves += 500;
534
535                        cout << "leaves=" << mVspStats.Leaves() << endl;
536                        Debug << "needed "
537                                  << TimeDiff(interTime, GetTime())*1e-3
538                                  << " secs to create 500 view cells" << endl;
539                        interTime = GetTime();
540                }
541
542                if (mCreatedViewCells == nViewCells)
543                {
544                        nViewCells += 500;
545
546                        cout << "generated " << mCreatedViewCells << " viewcells" << endl;
547                }
548        }
549
550        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
551        cout << "finished\n";
552
553        mVspStats.Stop();
554}
555
556
557bool VspOspTree::LocalTerminationCriteriaMet(const VspOspTraversalData &data) const
558{
559        return
560                (((int)data.mRays->size() <= mTermMinRays) ||
561                 (data.mPvs <= mTermMinPvs)   ||
562                 (data.mProbability <= mTermMinProbability) ||
563                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
564                 (data.mDepth >= mTermMaxDepth));
565}
566
567
568bool VspOspTree::GlobalTerminationCriteriaMet(const VspOspTraversalData &data) const
569{
570        return
571                (mOutOfMemory ||
572                (mVspStats.Leaves() >= mMaxViewCells) ||
573                (mGlobalCostMisses >= mTermGlobalCostMissTolerance));
574}
575
576
577// subdivide using a split plane queue
578VspNode *VspOspTree::Subdivide(VspOspSplitQueue &tQueue,
579                                                           VspOspSplitCandidate &splitCandidate)
580{
581        VspOspTraversalData &tData = splitCandidate.mParentData;
582
583        VspNode *newNode = tData.mNode;
584
585        if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
586        {       
587                VspOspTraversalData tFrontData;
588                VspOspTraversalData tBackData;
589
590                //-- continue subdivision
591               
592                // create new interior node and two leaf node
593                const AxisAlignedPlane splitPlane = splitCandidate.mSplitPlane;
594
595                newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData);
596       
597                const int maxCostMisses = splitCandidate.mMaxCostMisses;
598
599
600                // how often was max cost ratio missed in this branch?
601                tFrontData.mMaxCostMisses = maxCostMisses;
602                tBackData.mMaxCostMisses = maxCostMisses;
603                       
604               
605                if (1)
606                {
607                        float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
608                        float cBack = (float)tBackData.mPvs * tBackData.mProbability;
609                        float cData = (float)tData.mPvs * tData.mProbability;
610
611                       
612                        float costDecr =
613                                (cFront + cBack - cData) / mBox.GetVolume();
614
615                        mTotalCost += costDecr;
616                        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
617
618                        mSubdivisionStats
619                                        << "#ViewCells\n" << mVspStats.Leaves() << endl
620                                        << "#RenderCostDecrease\n" << -costDecr << endl
621                                        << "#SplitCandidateCost\n" << splitCandidate.GetCost() << endl
622                                        << "#TotalRenderCost\n" << mTotalCost << endl
623                                        << "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mVspStats.Leaves() << endl;
624                }
625
626       
627                //-- push the new split candidates on the stack
628                VspOspSplitCandidate frontCandidate;
629                VspOspSplitCandidate backCandidate;
630
631                EvalSplitCandidate(tFrontData, frontCandidate);
632                EvalSplitCandidate(tBackData, backCandidate);
633       
634                tQueue.push(frontCandidate);
635                tQueue.push(backCandidate);
636       
637                // delete old leaf node
638                DEL_PTR(tData.mNode);
639        }
640
641
642        //-- terminate traversal and create new view cell
643        if (newNode->IsLeaf())
644        {
645                VspLeaf *leaf = dynamic_cast<VspLeaf *>(newNode);
646                VspViewCell *viewCell = new VspViewCell();
647
648        leaf->SetViewCell(viewCell);
649               
650                //-- update pvs
651                int conSamp = 0;
652                float sampCon = 0.0f;
653                AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
654
655                // update scalar pvs size value
656                mViewCellsManager->SetScalarPvsSize(viewCell, viewCell->GetPvs().GetSize());
657
658                mVspStats.contributingSamples += conSamp;
659                mVspStats.sampleContributions +=(int) sampCon;
660
661                //-- store additional info
662                if (mStoreRays)
663                {
664                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
665                        for (it = tData.mRays->begin(); it != it_end; ++ it)
666                        {
667                                (*it).mRay->Ref();                     
668                                leaf->mVssRays.push_back((*it).mRay);
669                        }
670                }
671
672                // should I check here?
673                viewCell->mLeaf = leaf;
674
675                viewCell->SetVolume(tData.mProbability);
676        leaf->mProbability = tData.mProbability;
677
678                EvaluateLeafStats(tData);               
679        }
680
681        //-- cleanup
682        tData.Clear();
683
684        return newNode;
685}
686
687
688void VspOspTree::EvalSplitCandidate(VspOspTraversalData &tData,
689                                                                        VspOspSplitCandidate &splitData)
690{
691        VspOspTraversalData frontData;
692        VspOspTraversalData backData;
693
694        VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
695
696        // compute locally best split plane
697        const bool success =
698                SelectPlane(tData, splitData.mSplitPlane,
699                                        frontData.mProbability, backData.mProbability);
700
701        //TODO
702        // compute global decrease in render cost
703        splitData.mRenderCost = EvalRenderCostDecrease(splitData.mSplitPlane, tData);
704        splitData.mParentData = tData;
705        splitData.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
706}
707
708
709VspInterior *VspOspTree::SubdivideNode(const AxisAlignedPlane &splitPlane,
710                                                                           VspOspTraversalData &tData,
711                                                                           VspOspTraversalData &frontData,
712                                                                           VspOspTraversalData &backData)
713{
714        VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
715       
716        //-- the front and back traversal data is filled with the new values
717        frontData.mDepth = tData.mDepth + 1;
718        frontData.mRays = new RayInfoContainer();
719       
720        backData.mDepth = tData.mDepth + 1;
721        backData.mRays = new RayInfoContainer();
722       
723        //-- subdivide rays
724        SplitRays(splitPlane,
725                          *tData.mRays,
726                          *frontData.mRays,
727                          *backData.mRays);
728
729        //-- compute pvs
730        frontData.mPvs = ComputePvsSize(*frontData.mRays);
731        backData.mPvs = ComputePvsSize(*backData.mRays);
732
733        // split front and back node geometry and compute area
734        tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition,
735                                                         frontData.mBoundingBox, backData.mBoundingBox);
736
737
738        frontData.mProbability = frontData.mBoundingBox.GetVolume();
739        backData.mProbability = tData.mProbability - frontData.mProbability;
740
741       
742    ///////////////////////////////////////////
743        // subdivide further
744
745        // store maximal and minimal depth
746        if (tData.mDepth > mVspStats.maxDepth)
747        {
748                Debug << "max depth increases to " << tData.mDepth << " at " << mVspStats.Leaves() << " leaves" << endl;
749                mVspStats.maxDepth = tData.mDepth;
750        }
751
752        mVspStats.nodes += 2;
753
754   
755        VspInterior *interior = new VspInterior(splitPlane);
756
757#ifdef _DEBUG
758        Debug << interior << endl;
759#endif
760
761
762        //-- create front and back leaf
763
764        VspInterior *parent = leaf->GetParent();
765
766        // replace a link from node's parent
767        if (parent)
768        {
769                parent->ReplaceChildLink(leaf, interior);
770                interior->SetParent(parent);
771        }
772        else // new root
773        {
774                mRoot = interior;
775        }
776
777        // and setup child links
778        interior->SetupChildLinks(new VspLeaf(interior), new VspLeaf(interior));
779
780        frontData.mNode = interior->GetFront();
781        backData.mNode = interior->GetBack();
782
783        interior->mTimeStamp = mTimeStamp ++;
784       
785        return interior;
786}
787
788
789void VspOspTree::AddToPvs(VspLeaf *leaf,
790                                                  const RayInfoContainer &rays,
791                                                  float &sampleContributions,
792                                                  int &contributingSamples)
793{
794        sampleContributions = 0;
795        contributingSamples = 0;
796 
797        RayInfoContainer::const_iterator it, it_end = rays.end();
798 
799        ViewCellLeaf *vc = leaf->GetViewCell();
800 
801        // add contributions from samples to the PVS
802        for (it = rays.begin(); it != it_end; ++ it)
803        {
804                float sc = 0.0f;
805                VssRay *ray = (*it).mRay;
806
807                bool madeContrib = false;
808                float contribution;
809
810                if (ray->mTerminationObject)
811                {
812                        if (vc->AddPvsSample(ray->mTerminationObject, ray->mPdf, contribution))
813                                madeContrib = true;
814                        sc += contribution;
815                }
816         
817                if (ray->mOriginObject)
818                {
819                        if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
820                                madeContrib = true;
821                        sc += contribution;
822                }
823         
824          sampleContributions += sc;
825
826          if (madeContrib)
827                  ++ contributingSamples;
828               
829          //leaf->mVssRays.push_back(ray);
830        }
831}
832
833
834void VspOspTree::SortSplitCandidates(const RayInfoContainer &rays,
835                                                                         const int axis,
836                                                                         float minBand,
837                                                                         float maxBand)
838{
839        mSplitCandidates->clear();
840
841        int requestedSize = 2 * (int)(rays.size());
842        // creates a sorted split candidates array
843        if (mSplitCandidates->capacity() > 500000 &&
844                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
845        {
846        delete mSplitCandidates;
847                mSplitCandidates = new vector<SortableEntry>;
848        }
849
850        mSplitCandidates->reserve(requestedSize);
851
852        float pos;
853
854        // float values => don't compare with exact values
855        if (0)
856        {
857                minBand += Limits::Small;
858                maxBand -= Limits::Small;
859        }
860
861        // insert all queries
862        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
863        {
864                const bool positive = (*ri).mRay->HasPosDir(axis);
865                               
866                pos = (*ri).ExtrapOrigin(axis);
867                // clamp to min / max band
868                if (0) ClipValue(pos, minBand, maxBand);
869               
870                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
871                                                                        pos, (*ri).mRay));
872
873                pos = (*ri).ExtrapTermination(axis);
874                // clamp to min / max band
875                if (0) ClipValue(pos, minBand, maxBand);
876
877                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
878                                                                        pos, (*ri).mRay));
879        }
880
881        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
882}
883
884
885float VspOspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
886                                                                                  const AxisAlignedBox3 &box,
887                                                                                  const int pvsSize,
888                                                                                  const int axis,
889                                          float &position)
890{
891        const float minBox = box.Min(axis);
892        const float maxBox = box.Max(axis);
893
894        const float sizeBox = maxBox - minBox;
895
896        const float minBand = minBox + mMinBand * sizeBox;
897        const float maxBand = minBox + mMaxBand * sizeBox;
898
899        SortSplitCandidates(rays, axis, minBand, maxBand);
900
901        // go through the lists, count the number of objects left and right
902        // and evaluate the following cost funcion:
903        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
904
905        int pvsl = 0;
906        int pvsr = pvsSize;
907
908        int pvsBack = pvsl;
909        int pvsFront = pvsr;
910
911        float sum = (float)pvsSize * sizeBox;
912        float minSum = 1e20f;
913
914       
915        // if no border can be found, take mid split
916        position = minBox + 0.5f * sizeBox;
917       
918        // the relative cost ratio
919        float ratio = /*Limits::Infinity;*/99999999.0f;
920        bool splitPlaneFound = false;
921
922        Intersectable::NewMail();
923
924        RayInfoContainer::const_iterator ri, ri_end = rays.end();
925
926        // set all object as belonging to the front pvs
927        for(ri = rays.begin(); ri != ri_end; ++ ri)
928        {
929                Intersectable *oObject = (*ri).mRay->mOriginObject;
930                Intersectable *tObject = (*ri).mRay->mTerminationObject;
931
932                if (oObject)
933                {
934                        if (!oObject->Mailed())
935                        {
936                                oObject->Mail();
937                                oObject->mCounter = 1;
938                        }
939                        else
940                        {
941                                ++ oObject->mCounter;
942                        }
943                }
944                if (tObject)
945                {
946                        if (!tObject->Mailed())
947                        {
948                                tObject->Mail();
949                                tObject->mCounter = 1;
950                        }
951                        else
952                        {
953                                ++ tObject->mCounter;
954                        }
955                }
956        }
957
958        Intersectable::NewMail();
959
960        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
961
962        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
963        {
964                VssRay *ray;
965                ray = (*ci).ray;
966               
967                Intersectable *oObject = ray->mOriginObject;
968                Intersectable *tObject = ray->mTerminationObject;
969               
970
971                switch ((*ci).type)
972                {
973                        case SortableEntry::ERayMin:
974                                {
975                                        if (oObject && !oObject->Mailed())
976                                        {
977                                                oObject->Mail();
978                                                ++ pvsl;
979                                        }
980                                        if (tObject && !tObject->Mailed())
981                                        {
982                                                tObject->Mail();
983                                                ++ pvsl;
984                                        }
985                                        break;
986                                }
987                        case SortableEntry::ERayMax:
988                                {
989                                        if (oObject)
990                                        {
991                                                if (-- oObject->mCounter == 0)
992                                                        -- pvsr;
993                                        }
994
995                                        if (tObject)
996                                        {
997                                                if (-- tObject->mCounter == 0)
998                                                        -- pvsr;
999                                        }
1000
1001                                        break;
1002                                }
1003                }
1004
1005               
1006               
1007                // Note: sufficient to compare size of bounding boxes of front and back side?
1008                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
1009                {
1010                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
1011
1012                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << endl;
1013                        //Debug << "cost= " << sum << endl;
1014
1015                        if (sum < minSum)
1016                        {
1017                                splitPlaneFound = true;
1018
1019                                minSum = sum;
1020                                position = (*ci).value;
1021                               
1022                                pvsBack = pvsl;
1023                                pvsFront = pvsr;
1024                        }
1025                }
1026        }
1027       
1028       
1029        // -- compute cost
1030        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1031        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1032
1033        const float pOverall = sizeBox;
1034
1035        const float pBack = position - minBox;
1036        const float pFront = maxBox - position;
1037       
1038        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
1039    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1040        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1041       
1042        const float oldRenderCost = penaltyOld * pOverall;
1043        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1044
1045        if (splitPlaneFound)
1046        {
1047                ratio = newRenderCost / (oldRenderCost + Limits::Small);
1048        }
1049        //if (axis != 1)
1050        //Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1051         //    <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
1052
1053        return ratio;
1054}
1055
1056
1057float VspOspTree::SelectPlane(const VspOspTraversalData &tData,
1058                                                          AxisAlignedPlane &plane,
1059                                                          float &pFront,
1060                                                          float &pBack)
1061{
1062        float nPosition[3];
1063        float nCostRatio[3];
1064        float nProbFront[3];
1065        float nProbBack[3];
1066
1067        // create bounding box of node geometry
1068        AxisAlignedBox3 box = tData.mBoundingBox;
1069               
1070        int sAxis = 0;
1071        int bestAxis = -1;
1072
1073
1074        // if we use some kind of specialised fixed axis
1075    const bool useSpecialAxis =
1076                mOnlyDrivingAxis || mUseRandomAxis || mCirculatingAxis;
1077
1078        if (mUseRandomAxis)
1079                sAxis = Random(3);
1080        else if (mCirculatingAxis)
1081                sAxis = (tData.mAxis + 1) % 3;
1082        else if (mOnlyDrivingAxis)
1083                sAxis = box.Size().DrivingAxis();
1084
1085
1086        for (int axis = 0; axis < 3 ; ++ axis)
1087        {
1088                if (!useSpecialAxis || (axis == sAxis))
1089                {
1090                        //-- place split plane using heuristics
1091
1092                        if (mUseCostHeuristics)
1093                        {
1094                                nCostRatio[axis] =
1095                                        BestCostRatioHeuristics(*tData.mRays,
1096                                                                                    box,
1097                                                                                        tData.mPvs,
1098                                                                                        axis,
1099                                                                                        nPosition[axis]);                       
1100                        }
1101                        else //-- split plane position is spatial median
1102                        {
1103                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1104
1105                                nCostRatio[axis] = EvalSplitCost(tData,
1106                                                                                                 box,
1107                                                                                                 axis,
1108                                                                                                 nPosition[axis],
1109                                                                                                 nProbFront[axis],
1110                                                                                                 nProbBack[axis]);
1111                        }
1112                                               
1113                        if (bestAxis == -1)
1114                        {
1115                                bestAxis = axis;
1116                        }
1117                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1118                        {
1119                                bestAxis = axis;
1120                        }
1121                }
1122        }
1123
1124        //-- assign values
1125        plane.mAxis = bestAxis;
1126        pFront = nProbFront[bestAxis];
1127        pBack = nProbBack[bestAxis];
1128
1129        //-- split plane position
1130    plane.mPosition = nPosition[bestAxis];
1131
1132        return nCostRatio[bestAxis];
1133}
1134
1135
1136inline void VspOspTree::GenerateUniqueIdsForPvs()
1137{
1138        Intersectable::NewMail(); sBackId = Intersectable::sMailId;
1139        Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
1140        Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
1141}
1142
1143
1144float VspOspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
1145                                                                                 const VspOspTraversalData &data) const
1146{
1147        float pvsFront = 0;
1148        float pvsBack = 0;
1149        float totalPvs = 0;
1150
1151        // probability that view point lies in back / front node
1152        float pOverall = data.mProbability;
1153        float pFront = 0;
1154        float pBack = 0;
1155
1156
1157        // create unique ids for pvs heuristics
1158        GenerateUniqueIdsForPvs();
1159       
1160        for (int i = 0; i < data.mRays->size(); ++ i)
1161        {
1162                RayInfo rayInf = (*data.mRays)[i];
1163
1164                float t;
1165                VssRay *ray = rayInf.mRay;
1166                const int cf =
1167                        rayInf.ComputeRayIntersection(candidatePlane.mAxis,
1168                                                                                  candidatePlane.mPosition, t);
1169
1170                // find front and back pvs for origing and termination object
1171                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1172                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1173        }
1174
1175        AxisAlignedBox3 frontBox;
1176        AxisAlignedBox3 backBox;
1177
1178        data.mBoundingBox.Split(candidatePlane.mAxis, candidatePlane.mPosition, frontBox, backBox);
1179
1180        pFront = frontBox.GetVolume();
1181        pBack = pOverall - pFront;
1182               
1183
1184        // -- pvs rendering heuristics
1185        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1186        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1187
1188        // only render cost heuristics or combined with standard deviation
1189        const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
1190    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1191        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1192                       
1193        const float oldRenderCost = pOverall * penaltyOld;
1194        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1195
1196        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
1197        const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBox.GetVolume();
1198       
1199
1200#if 1
1201        // take render cost of node into account to avoid being stuck in a local minimum
1202        const float normalizedOldRenderCost = oldRenderCost / mBox.GetVolume();
1203        return 0.99f * renderCostDecrease + 0.01f * normalizedOldRenderCost;
1204#else
1205        return renderCostDecrease;
1206#endif
1207}
1208
1209
1210
1211float VspOspTree::EvalSplitCost(const VspOspTraversalData &data,
1212                                                                const AxisAlignedBox3 &box,
1213                                                                const int axis,
1214                                                                const float &position,                                                                           
1215                                                                float &pFront,
1216                                                                float &pBack) const
1217{
1218        float pvsTotal = 0;
1219        float pvsFront = 0;
1220        float pvsBack = 0;
1221       
1222        // create unique ids for pvs heuristics
1223        GenerateUniqueIdsForPvs();
1224
1225        const int pvsSize = data.mPvs;
1226
1227        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1228
1229        // this is the main ray classification loop!
1230        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
1231        {
1232                // determine the side of this ray with respect to the plane
1233                float t;
1234                const int side = (*rit).ComputeRayIntersection(axis, position, t);
1235       
1236                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
1237                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
1238        }
1239
1240        //-- pvs heuristics
1241        float pOverall;
1242
1243        //-- compute heurstics
1244        //-- we take simplified computation for mid split
1245               
1246        pOverall = data.mProbability;
1247        pBack = pFront = pOverall * 0.5f;
1248       
1249       
1250#ifdef _DEBUG
1251        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
1252        Debug << pFront << " " << pBack << " " << pOverall << endl;
1253#endif
1254
1255       
1256        const float newCost = pvsBack * pBack + pvsFront * pFront;
1257        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
1258
1259        return  (mCtDivCi + newCost) / oldCost;
1260}
1261
1262
1263void VspOspTree::AddObjToPvs(Intersectable *obj,
1264                                                         const int cf,
1265                                                         float &frontPvs,
1266                                                         float &backPvs,
1267                                                         float &totalPvs) const
1268{
1269        if (!obj)
1270                return;
1271
1272        const float renderCost = mViewCellsManager->EvalRenderCost(obj);
1273
1274        // new object
1275        if ((obj->mMailbox != sFrontId) &&
1276                (obj->mMailbox != sBackId) &&
1277                (obj->mMailbox != sFrontAndBackId))
1278        {
1279                totalPvs += renderCost;
1280        }
1281
1282        // TODO: does this really belong to no pvs?
1283        //if (cf == Ray::COINCIDENT) return;
1284
1285        // object belongs to both PVS
1286        if (cf >= 0)
1287        {
1288                if ((obj->mMailbox != sFrontId) &&
1289                        (obj->mMailbox != sFrontAndBackId))
1290                {
1291                        frontPvs += renderCost;
1292               
1293                        if (obj->mMailbox == sBackId)
1294                                obj->mMailbox = sFrontAndBackId;
1295                        else
1296                                obj->mMailbox = sFrontId;
1297                }
1298        }
1299
1300        if (cf <= 0)
1301        {
1302                if ((obj->mMailbox != sBackId) &&
1303                        (obj->mMailbox != sFrontAndBackId))
1304                {
1305                        backPvs += renderCost;
1306               
1307                        if (obj->mMailbox == sFrontId)
1308                                obj->mMailbox = sFrontAndBackId;
1309                        else
1310                                obj->mMailbox = sBackId;
1311                }
1312        }
1313}
1314
1315
1316void VspOspTree::CollectLeaves(vector<VspLeaf *> &leaves,
1317                                                           const bool onlyUnmailed,
1318                                                           const int maxPvsSize) const
1319{
1320        stack<VspNode *> nodeStack;
1321        nodeStack.push(mRoot);
1322
1323        while (!nodeStack.empty())
1324        {
1325                VspNode *node = nodeStack.top();
1326                nodeStack.pop();
1327               
1328                if (node->IsLeaf())
1329                {
1330                        // test if this leaf is in valid view space
1331                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1332                        if (leaf->TreeValid() &&
1333                                (!onlyUnmailed || !leaf->Mailed()) &&
1334                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
1335                        {
1336                                leaves.push_back(leaf);
1337                        }
1338                }
1339                else
1340                {
1341                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1342
1343                        nodeStack.push(interior->GetBack());
1344                        nodeStack.push(interior->GetFront());
1345                }
1346        }
1347}
1348
1349
1350AxisAlignedBox3 VspOspTree::GetBoundingBox() const
1351{
1352        return mBox;
1353}
1354
1355
1356VspNode *VspOspTree::GetRoot() const
1357{
1358        return mRoot;
1359}
1360
1361
1362void VspOspTree::EvaluateLeafStats(const VspOspTraversalData &data)
1363{
1364        // the node became a leaf -> evaluate stats for leafs
1365        VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);
1366
1367
1368        if (data.mPvs > mVspStats.maxPvs)
1369        {
1370                mVspStats.maxPvs = data.mPvs;
1371        }
1372
1373        mVspStats.pvs += data.mPvs;
1374
1375        if (data.mDepth < mVspStats.minDepth)
1376        {
1377                mVspStats.minDepth = data.mDepth;
1378        }
1379       
1380        if (data.mDepth >= mTermMaxDepth)
1381        {
1382        ++ mVspStats.maxDepthNodes;
1383                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
1384        }
1385
1386        // accumulate rays to compute rays /  leaf
1387        mVspStats.accumRays += (int)data.mRays->size();
1388
1389        if (data.mPvs < mTermMinPvs)
1390                ++ mVspStats.minPvsNodes;
1391
1392        if ((int)data.mRays->size() < mTermMinRays)
1393                ++ mVspStats.minRaysNodes;
1394
1395        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1396                ++ mVspStats.maxRayContribNodes;
1397
1398        if (data.mProbability <= mTermMinProbability)
1399                ++ mVspStats.minProbabilityNodes;
1400       
1401        // accumulate depth to compute average depth
1402        mVspStats.accumDepth += data.mDepth;
1403
1404        ++ mCreatedViewCells;
1405
1406#ifdef _DEBUG
1407        Debug << "BSP stats: "
1408                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1409                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1410          //              << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
1411                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1412                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
1413                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1414#endif
1415}
1416
1417
1418void VspOspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
1419{
1420        ViewCell::NewMail();
1421        CollectViewCells(mRoot, onlyValid, viewCells, true);
1422}
1423
1424
1425void VspOspTree::CollapseViewCells()
1426{
1427// TODO
1428#if HAS_TO_BE_REDONE
1429        stack<VspNode *> nodeStack;
1430
1431        if (!mRoot)
1432                return;
1433
1434        nodeStack.push(mRoot);
1435       
1436        while (!nodeStack.empty())
1437        {
1438                VspNode *node = nodeStack.top();
1439                nodeStack.pop();
1440               
1441                if (node->IsLeaf())
1442        {
1443                        BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1444
1445                        if (!viewCell->GetValid())
1446                        {
1447                                BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1448       
1449                                ViewCellContainer leaves;
1450                                mViewCellsTree->CollectLeaves(viewCell, leaves);
1451
1452                                ViewCellContainer::const_iterator it, it_end = leaves.end();
1453
1454                                for (it = leaves.begin(); it != it_end; ++ it)
1455                                {
1456                                        VspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
1457                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
1458                                        ++ mVspStats.invalidLeaves;
1459                                }
1460
1461                                // add to unbounded view cell
1462                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
1463                                DEL_PTR(viewCell);
1464                        }
1465                }
1466                else
1467                {
1468                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1469               
1470                        nodeStack.push(interior->GetFront());
1471                        nodeStack.push(interior->GetBack());
1472                }
1473        }
1474
1475        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1476#endif
1477}
1478
1479
1480void VspOspTree::CollectRays(VssRayContainer &rays)
1481{
1482        vector<VspLeaf *> leaves;
1483
1484        vector<VspLeaf *>::const_iterator lit, lit_end = leaves.end();
1485
1486        for (lit = leaves.begin(); lit != lit_end; ++ lit)
1487        {
1488                VspLeaf *leaf = *lit;
1489                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
1490
1491                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
1492                        rays.push_back(*rit);
1493        }
1494}
1495
1496
1497void VspOspTree::ValidateTree()
1498{
1499        stack<VspNode *> nodeStack;
1500
1501        if (!mRoot)
1502                return;
1503
1504        nodeStack.push(mRoot);
1505       
1506        mVspStats.invalidLeaves = 0;
1507        while (!nodeStack.empty())
1508        {
1509                VspNode *node = nodeStack.top();
1510                nodeStack.pop();
1511               
1512                if (node->IsLeaf())
1513                {
1514                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1515
1516                        if (!leaf->GetViewCell()->GetValid())
1517                                ++ mVspStats.invalidLeaves;
1518
1519                        // validity flags don't match => repair
1520                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
1521                        {
1522                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
1523                                PropagateUpValidity(leaf);
1524                        }
1525                }
1526                else
1527                {
1528                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1529               
1530                        nodeStack.push(interior->GetFront());
1531                        nodeStack.push(interior->GetBack());
1532                }
1533        }
1534
1535        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1536}
1537
1538
1539
1540void VspOspTree::CollectViewCells(VspNode *root,
1541                                                                  bool onlyValid,
1542                                                                  ViewCellContainer &viewCells,
1543                                                                  bool onlyUnmailed) const
1544{
1545        stack<VspNode *> nodeStack;
1546
1547        if (!root)
1548                return;
1549
1550        nodeStack.push(root);
1551       
1552        while (!nodeStack.empty())
1553        {
1554                VspNode *node = nodeStack.top();
1555                nodeStack.pop();
1556               
1557                if (node->IsLeaf())
1558                {
1559                        if (!onlyValid || node->TreeValid())
1560                        {
1561                                ViewCellLeaf *leafVc = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1562
1563                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
1564                                               
1565                                if (!onlyUnmailed || !viewCell->Mailed())
1566                                {
1567                                        viewCell->Mail();
1568                                        viewCells.push_back(viewCell);
1569                                }
1570                        }
1571                }
1572                else
1573                {
1574                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1575               
1576                        nodeStack.push(interior->GetFront());
1577                        nodeStack.push(interior->GetBack());
1578                }
1579        }
1580
1581}
1582
1583
1584int VspOspTree::FindNeighbors(VspLeaf *n,
1585                                                          vector<VspLeaf *> &neighbors,
1586                                                          const bool onlyUnmailed) const
1587{
1588        stack<VspNode *> nodeStack;
1589
1590        nodeStack.push(mRoot);
1591
1592        const AxisAlignedBox3 box = GetBBox(n);
1593
1594        while (!nodeStack.empty())
1595        {
1596                VspNode *node = nodeStack.top();
1597                nodeStack.pop();
1598
1599                if (node->IsLeaf())
1600                {
1601                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1602
1603                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
1604                                neighbors.push_back(leaf);
1605                }
1606                else
1607                {
1608                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1609                       
1610                        if (interior->GetPosition() > box.Max(interior->GetAxis()))
1611                                nodeStack.push(interior->GetBack());
1612                        else
1613                        {
1614                                if (interior->GetPosition() < box.Min(interior->GetAxis()))
1615                                        nodeStack.push(interior->GetFront());
1616                                else
1617                                {
1618                                        // random decision
1619                                        nodeStack.push(interior->GetBack());
1620                                        nodeStack.push(interior->GetFront());
1621                                }
1622                        }
1623                }
1624        }
1625
1626        return (int)neighbors.size();
1627}
1628
1629
1630// Find random neighbor which was not mailed
1631VspLeaf *VspOspTree::GetRandomLeaf(const Plane3 &plane)
1632{
1633        stack<VspNode *> nodeStack;
1634        nodeStack.push(mRoot);
1635 
1636        int mask = rand();
1637 
1638        while (!nodeStack.empty())
1639        {
1640                VspNode *node = nodeStack.top();
1641               
1642                nodeStack.pop();
1643               
1644                if (node->IsLeaf())
1645                {
1646                        return dynamic_cast<VspLeaf *>(node);
1647                }
1648                else
1649                {
1650                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1651                        VspNode *next;
1652                       
1653                        if (GetBBox(interior->GetBack()).Side(plane) < 0)
1654                                next = interior->GetFront();
1655            else
1656                                if (GetBBox(interior->GetFront()).Side(plane) < 0)
1657                                        next = interior->GetBack();
1658                                else
1659                                {
1660                                        // random decision
1661                                        if (mask & 1)
1662                                                next = interior->GetBack();
1663                                        else
1664                                                next = interior->GetFront();
1665                                        mask = mask >> 1;
1666                                }
1667
1668                                nodeStack.push(next);
1669                }
1670        }
1671 
1672        return NULL;
1673}
1674
1675
1676VspLeaf *VspOspTree::GetRandomLeaf(const bool onlyUnmailed)
1677{
1678        stack<VspNode *> nodeStack;
1679
1680        nodeStack.push(mRoot);
1681
1682        int mask = rand();
1683
1684        while (!nodeStack.empty())
1685        {
1686                VspNode *node = nodeStack.top();
1687                nodeStack.pop();
1688
1689                if (node->IsLeaf())
1690                {
1691                        if ( (!onlyUnmailed || !node->Mailed()) )
1692                                return dynamic_cast<VspLeaf *>(node);
1693                }
1694                else
1695                {
1696                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1697
1698                        // random decision
1699                        if (mask & 1)
1700                                nodeStack.push(interior->GetBack());
1701                        else
1702                                nodeStack.push(interior->GetFront());
1703
1704                        mask = mask >> 1;
1705                }
1706        }
1707
1708        return NULL;
1709}
1710
1711
1712int VspOspTree::ComputePvsSize(const RayInfoContainer &rays) const
1713{
1714        int pvsSize = 0;
1715
1716        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1717
1718        Intersectable::NewMail();
1719
1720        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1721        {
1722                VssRay *ray = (*rit).mRay;
1723
1724                if (ray->mOriginObject)
1725                {
1726                        if (!ray->mOriginObject->Mailed())
1727                        {
1728                                ray->mOriginObject->Mail();
1729                                ++ pvsSize;
1730                        }
1731                }
1732                if (ray->mTerminationObject)
1733                {
1734                        if (!ray->mTerminationObject->Mailed())
1735                        {
1736                                ray->mTerminationObject->Mail();
1737                                ++ pvsSize;
1738                        }
1739                }
1740        }
1741
1742        return pvsSize;
1743}
1744
1745
1746float VspOspTree::GetEpsilon() const
1747{
1748        return mEpsilon;
1749}
1750
1751
1752int VspOspTree::CastLineSegment(const Vector3 &origin,
1753                                                                const Vector3 &termination,
1754                                                                ViewCellContainer &viewcells)
1755{
1756        int hits = 0;
1757
1758        float mint = 0.0f, maxt = 1.0f;
1759        const Vector3 dir = termination - origin;
1760
1761        stack<LineTraversalData> tStack;
1762
1763        Intersectable::NewMail();
1764        ViewCell::NewMail();
1765
1766        Vector3 entp = origin;
1767        Vector3 extp = termination;
1768
1769        VspNode *node = mRoot;
1770        VspNode *farChild;
1771
1772        float position;
1773        int axis;
1774
1775        while (1)
1776        {
1777                if (!node->IsLeaf())
1778                {
1779                        VspInterior *in = dynamic_cast<VspInterior *>(node);
1780                        position = in->GetPosition();
1781                        axis = in->GetAxis();
1782
1783                        if (entp[axis] <= position)
1784                        {
1785                                if (extp[axis] <= position)
1786                                {
1787                                        node = in->GetBack();
1788                                        // cases N1,N2,N3,P5,Z2,Z3
1789                                        continue;
1790                                } else
1791                                {
1792                                        // case N4
1793                                        node = in->GetBack();
1794                                        farChild = in->GetFront();
1795                                }
1796                        }
1797                        else
1798                        {
1799                                if (position <= extp[axis])
1800                                {
1801                                        node = in->GetFront();
1802                                        // cases P1,P2,P3,N5,Z1
1803                                        continue;
1804                                }
1805                                else
1806                                {
1807                                        node = in->GetFront();
1808                                        farChild = in->GetBack();
1809                                        // case P4
1810                                }
1811                        }
1812
1813                        // $$ modification 3.5.2004 - hints from Kamil Ghais
1814                        // case N4 or P4
1815                        const float tdist = (position - origin[axis]) / dir[axis];
1816                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
1817
1818                        extp = origin + dir * tdist;
1819                        maxt = tdist;
1820                }
1821                else
1822                {
1823                        // compute intersection with all objects in this leaf
1824                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1825                        ViewCell *vc = leaf->GetViewCell();
1826
1827                        if (!vc->Mailed())
1828                        {
1829                                vc->Mail();
1830                                viewcells.push_back(vc);
1831                                ++ hits;
1832                        }
1833#if 0
1834                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
1835#endif
1836                        // get the next node from the stack
1837                        if (tStack.empty())
1838                                break;
1839
1840                        entp = extp;
1841                        mint = maxt;
1842                       
1843                        LineTraversalData &s  = tStack.top();
1844                        node = s.mNode;
1845                        extp = s.mExitPoint;
1846                        maxt = s.mMaxT;
1847                        tStack.pop();
1848                }
1849        }
1850
1851        return hits;
1852}
1853
1854
1855int VspOspTree::CastRay(Ray &ray)
1856{
1857        int hits = 0;
1858
1859        stack<LineTraversalData> tStack;
1860        const Vector3 dir = ray.GetDir();
1861
1862        float maxt, mint;
1863
1864        if (!mBox.GetRaySegment(ray, mint, maxt))
1865                return 0;
1866
1867        Intersectable::NewMail();
1868        ViewCell::NewMail();
1869
1870        Vector3 entp = ray.Extrap(mint);
1871        Vector3 extp = ray.Extrap(maxt);
1872
1873        const Vector3 origin = entp;
1874
1875        VspNode *node = mRoot;
1876        VspNode *farChild = NULL;
1877
1878        float position;
1879        int axis;
1880
1881        while (1)
1882        {
1883                if (!node->IsLeaf())
1884                {
1885                        VspInterior *in = dynamic_cast<VspInterior *>(node);
1886                        position = in->GetPosition();
1887                        axis = in->GetAxis();
1888
1889                        if (entp[axis] <= position)
1890                        {
1891                                if (extp[axis] <= position)
1892                                {
1893                                        node = in->GetBack();
1894                                        // cases N1,N2,N3,P5,Z2,Z3
1895                                        continue;
1896                                } else
1897                                {
1898                                        // case N4
1899                                        node = in->GetBack();
1900                                        farChild = in->GetFront();
1901                                }
1902                        }
1903                        else
1904                        {
1905                                if (position <= extp[axis])
1906                                {
1907                                        node = in->GetFront();
1908                                        // cases P1,P2,P3,N5,Z1
1909                                        continue;
1910                                }
1911                                else
1912                                {
1913                                        node = in->GetFront();
1914                                        farChild = in->GetBack();
1915                                        // case P4
1916                                }
1917                        }
1918
1919                        // $$ modification 3.5.2004 - hints from Kamil Ghais
1920                        // case N4 or P4
1921                        const float tdist = (position - origin[axis]) / dir[axis];
1922                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
1923                        extp = origin + dir * tdist;
1924                        maxt = tdist;
1925                }
1926                else
1927                {
1928                        // compute intersection with all objects in this leaf
1929                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1930                        ViewCell *vc = leaf->GetViewCell();
1931
1932                        if (!vc->Mailed())
1933                        {
1934                                vc->Mail();
1935                                // todo: add view cells to ray
1936                                ++ hits;
1937                        }
1938
1939#if 0
1940                                leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
1941#endif
1942                        // get the next node from the stack
1943                        if (tStack.empty())
1944                                break;
1945
1946                        entp = extp;
1947                        mint = maxt;
1948                       
1949                        LineTraversalData &s  = tStack.top();
1950                        node = s.mNode;
1951                        extp = s.mExitPoint;
1952                        maxt = s.mMaxT;
1953                        tStack.pop();
1954                }
1955        }
1956
1957        return hits;
1958}
1959
1960
1961VspNode *VspOspTree::CollapseTree(VspNode *node, int &collapsed)
1962{
1963// TODO
1964#if HAS_TO_BE_REDONE
1965        if (node->IsLeaf())
1966                return node;
1967
1968        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1969
1970        VspNode *front = CollapseTree(interior->GetFront(), collapsed);
1971        VspNode *back = CollapseTree(interior->GetBack(), collapsed);
1972
1973        if (front->IsLeaf() && back->IsLeaf())
1974        {
1975                VspLeaf *frontLeaf = dynamic_cast<VspLeaf *>(front);
1976                VspLeaf *backLeaf = dynamic_cast<VspLeaf *>(back);
1977
1978                //-- collapse tree
1979                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
1980                {
1981                        BspViewCell *vc = frontLeaf->GetViewCell();
1982
1983                        VspLeaf *leaf = new VspLeaf(interior->GetParent(), vc);
1984                        leaf->SetTreeValid(frontLeaf->TreeValid());
1985
1986                        // replace a link from node's parent
1987                        if (leaf->GetParent())
1988                                leaf->GetParent()->ReplaceChildLink(node, leaf);
1989                        else
1990                                mRoot = leaf;
1991
1992                        ++ collapsed;
1993                        delete interior;
1994
1995                        return leaf;
1996                }
1997        }
1998#endif
1999        return node;
2000}
2001
2002
2003int VspOspTree::CollapseTree()
2004{
2005        int collapsed = 0;
2006        //TODO
2007#if HAS_TO_BE_REDONE
2008        (void)CollapseTree(mRoot, collapsed);
2009
2010        // revalidate leaves
2011        RepairViewCellsLeafLists();
2012#endif
2013        return collapsed;
2014}
2015
2016
2017void VspOspTree::RepairViewCellsLeafLists()
2018{
2019// TODO
2020#if HAS_TO_BE_REDONE
2021        // list not valid anymore => clear
2022        stack<VspNode *> nodeStack;
2023        nodeStack.push(mRoot);
2024
2025        ViewCell::NewMail();
2026
2027        while (!nodeStack.empty())
2028        {
2029                VspNode *node = nodeStack.top();
2030                nodeStack.pop();
2031
2032                if (node->IsLeaf())
2033                {
2034                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2035
2036                        BspViewCell *viewCell = leaf->GetViewCell();
2037
2038                        if (!viewCell->Mailed())
2039                        {
2040                                viewCell->mLeaves.clear();
2041                                viewCell->Mail();
2042                        }
2043       
2044                        viewCell->mLeaves.push_back(leaf);
2045
2046                }
2047                else
2048                {
2049                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2050
2051                        nodeStack.push(interior->GetFront());
2052                        nodeStack.push(interior->GetBack());
2053                }
2054        }
2055// TODO
2056#endif
2057}
2058
2059
2060
2061ViewCell *VspOspTree::GetViewCell(const Vector3 &point, const bool active)
2062{
2063        if (mRoot == NULL)
2064                return NULL;
2065
2066        stack<VspNode *> nodeStack;
2067        nodeStack.push(mRoot);
2068 
2069        ViewCellLeaf *viewcell = NULL;
2070 
2071        while (!nodeStack.empty()) 
2072        {
2073                VspNode *node = nodeStack.top();
2074                nodeStack.pop();
2075       
2076                if (node->IsLeaf())
2077                {
2078                        viewcell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
2079                        break;
2080                }
2081                else   
2082                {       
2083                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2084     
2085                        // random decision
2086                        if (interior->GetPosition() - point[interior->GetAxis()] < 0)
2087                                nodeStack.push(interior->GetBack());
2088                        else
2089                                nodeStack.push(interior->GetFront());
2090                }
2091        }
2092 
2093        if (active)
2094                return mViewCellsTree->GetActiveViewCell(viewcell);
2095        else
2096                return viewcell;
2097}
2098
2099
2100bool VspOspTree::ViewPointValid(const Vector3 &viewPoint) const
2101{
2102        VspNode *node = mRoot;
2103
2104        while (1)
2105        {
2106                // early exit
2107                if (node->TreeValid())
2108                        return true;
2109
2110                if (node->IsLeaf())
2111                        return false;
2112                       
2113                VspInterior *in = dynamic_cast<VspInterior *>(node);
2114                                       
2115                if (in->GetPosition() - viewPoint[in->GetAxis()] <= 0)
2116                {
2117                        node = in->GetBack();
2118                }
2119                else
2120                {
2121                        node = in->GetFront();
2122                }
2123        }
2124
2125        // should never come here
2126        return false;
2127}
2128
2129
2130void VspOspTree::PropagateUpValidity(VspNode *node)
2131{
2132        const bool isValid = node->TreeValid();
2133
2134        // propagative up invalid flag until only invalid nodes exist over this node
2135        if (!isValid)
2136        {
2137                while (!node->IsRoot() && node->GetParent()->TreeValid())
2138                {
2139                        node = node->GetParent();
2140                        node->SetTreeValid(false);
2141                }
2142        }
2143        else
2144        {
2145                // propagative up valid flag until one of the subtrees is invalid
2146                while (!node->IsRoot() && !node->TreeValid())
2147                {
2148            node = node->GetParent();
2149                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2150                       
2151                        // the parent is valid iff both leaves are valid
2152                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
2153                                                           interior->GetFront()->TreeValid());
2154                }
2155        }
2156}
2157
2158#if ZIPPED_VIEWCELLS
2159bool VspOspTree::Export(ogzstream &stream)
2160#else
2161bool VspOspTree::Export(ofstream &stream)
2162#endif
2163{
2164        ExportNode(mRoot, stream);
2165
2166        return true;
2167}
2168
2169
2170#if ZIPPED_VIEWCELLS
2171void VspOspTree::ExportNode(VspNode *node, ogzstream &stream)
2172#else
2173void VspOspTree::ExportNode(VspNode *node, ofstream &stream)
2174#endif
2175{
2176        if (node->IsLeaf())
2177        {
2178                VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2179                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
2180
2181                int id = -1;
2182                if (viewCell != mOutOfBoundsCell)
2183                        id = viewCell->GetId();
2184
2185                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
2186        }
2187        else
2188        {       
2189                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2190       
2191                AxisAlignedPlane plane = interior->GetPlane();
2192                stream << "<Interior plane=\"" << plane.mPosition << " "
2193                           << plane.mAxis << "\">" << endl;
2194
2195                ExportNode(interior->GetBack(), stream);
2196                ExportNode(interior->GetFront(), stream);
2197
2198                stream << "</Interior>" << endl;
2199        }
2200}
2201
2202
2203int VspOspTree::SplitRays(const AxisAlignedPlane &plane,
2204                                                  RayInfoContainer &rays,
2205                                                  RayInfoContainer &frontRays,
2206                                                  RayInfoContainer &backRays) const
2207{
2208        int splits = 0;
2209
2210        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2211
2212        for (rit = rays.begin(); rit != rit_end; ++ rit)
2213        {
2214                RayInfo bRay = *rit;
2215               
2216                VssRay *ray = bRay.mRay;
2217                float t;
2218
2219                // get classification and receive new t
2220                //-- test if start point behind or in front of plane
2221                const int side = plane.ComputeRayIntersection(bRay, t);
2222
2223                if (side == 0)
2224                {
2225                        ++ splits;
2226
2227                        if (ray->HasPosDir(plane.mAxis))
2228                        {
2229                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2230                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2231                        }
2232                        else
2233                        {
2234                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2235                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2236                        }
2237                }
2238                else if (side == 1)
2239                {
2240                        frontRays.push_back(bRay);
2241                }
2242                else
2243                {
2244                        backRays.push_back(bRay);
2245                }
2246
2247        }
2248
2249        return splits;
2250}
2251
2252
2253AxisAlignedBox3 VspOspTree::GetBBox(VspNode *node) const
2254{
2255        if (!node->GetParent())
2256                return mBox;
2257
2258        if (!node->IsLeaf())
2259        {
2260                return (dynamic_cast<VspInterior *>(node))->GetBox();           
2261        }
2262
2263        VspInterior *parent = dynamic_cast<VspInterior *>(node->GetParent());
2264
2265        AxisAlignedBox3 box(parent->GetBox());
2266
2267        if (parent->GetFront() == node)
2268                box.SetMin(parent->GetAxis(), parent->GetPosition());
2269        else
2270                box.SetMax(parent->GetAxis(), parent->GetPosition());
2271
2272        return box;
2273}
2274
2275
2276}
Note: See TracBrowser for help on using the repository browser.