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

Revision 1077, 78.7 KB checked in by mattausch, 18 years ago (diff)

worked on object space /view space partitioning
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#include "KdTree.h"
19
20
21namespace GtpVisibilityPreprocessor {
22
23#define USE_FIXEDPOINT_T 0
24
25
26//-- static members
27
28int VspTree::sFrontId = 0;
29int VspTree::sBackId = 0;
30int VspTree::sFrontAndBackId = 0;
31
32
33
34// pvs penalty can be different from pvs size
35inline static float EvalPvsPenalty(const int pvs,
36                                                                   const int lower,
37                                                                   const int upper)
38{
39        // clamp to minmax values
40        if (pvs < lower)
41                return (float)lower;
42        else if (pvs > upper)
43                return (float)upper;
44
45        return (float)pvs;
46}
47
48
49int VspNode::sMailId = 1;
50
51
52
53void VspTreeStatistics::Print(ostream &app) const
54{
55        app << "===== VspTree statistics ===============\n";
56
57        app << setprecision(4);
58
59        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
60
61        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
62
63        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
64
65        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
66
67        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits[0] + splits[1] + splits[2] << endl;
68
69        app << "#N_SPLITS ( Number of splits in axes x y z)\n";
70
71        for (int i = 0; i < 3; ++ i)
72                app << splits[i] << " ";
73        app << endl;
74
75        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
76                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
77
78        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n"
79                << minPvsNodes * 100 / (double)Leaves() << endl;
80
81        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n"
82                << minRaysNodes * 100 / (double)Leaves() << endl;
83
84        app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
85                << maxCostNodes * 100 / (double)Leaves() << endl;
86
87        app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
88                << minProbabilityNodes * 100 / (double)Leaves() << endl;
89
90        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"
91                <<      maxRayContribNodes * 100 / (double)Leaves() << endl;
92
93        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
94
95        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
96
97        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
98
99        app << "#N_INVALIDLEAVES (number of invalid leaves )\n" << invalidLeaves << endl;
100
101        app << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
102        //app << "#N_PVS: " << pvs << endl;
103
104        app << "===== END OF VspTree statistics ==========\n";
105}
106
107
108/******************************************************************/
109/*                  class VspNode implementation                  */
110/******************************************************************/
111
112
113VspNode::VspNode():
114mParent(NULL), mTreeValid(true), mTimeStamp(0)
115{}
116
117
118VspNode::VspNode(VspInterior *parent):
119mParent(parent), mTreeValid(true)
120{}
121
122
123bool VspNode::IsRoot() const
124{
125        return mParent == NULL;
126}
127
128
129VspInterior *VspNode::GetParent()
130{
131        return mParent;
132}
133
134
135void VspNode::SetParent(VspInterior *parent)
136{
137        mParent = parent;
138}
139
140
141bool VspNode::IsSibling(VspNode *n) const
142{
143        return  ((this != n) && mParent &&
144                         (mParent->GetFront() == n) || (mParent->GetBack() == n));
145}
146
147
148int VspNode::GetDepth() const
149{
150        int depth = 0;
151        VspNode *p = mParent;
152       
153        while (p)
154        {
155                p = p->mParent;
156                ++ depth;
157        }
158
159        return depth;
160}
161
162
163bool VspNode::TreeValid() const
164{
165        return mTreeValid;
166}
167
168
169void VspNode::SetTreeValid(const bool v)
170{
171        mTreeValid = v;
172}
173
174
175/****************************************************************/
176/*              class VspInterior implementation                */
177/****************************************************************/
178
179
180VspInterior::VspInterior(const AxisAlignedPlane &plane):
181mPlane(plane), mFront(NULL), mBack(NULL)
182{}
183
184
185VspInterior::~VspInterior()
186{
187        DEL_PTR(mFront);
188        DEL_PTR(mBack);
189}
190
191
192bool VspInterior::IsLeaf() const
193{
194        return false;
195}
196
197
198VspNode *VspInterior::GetBack()
199{
200        return mBack;
201}
202
203
204VspNode *VspInterior::GetFront()
205{
206        return mFront;
207}
208
209
210AxisAlignedPlane VspInterior::GetPlane() const
211{
212        return mPlane;
213}
214
215
216float VspInterior::GetPosition() const
217{
218        return mPlane.mPosition;
219}
220
221
222int VspInterior::GetAxis() const
223{
224        return mPlane.mAxis;
225}
226
227
228void VspInterior::ReplaceChildLink(VspNode *oldChild, VspNode *newChild)
229{
230        if (mBack == oldChild)
231                mBack = newChild;
232        else
233                mFront = newChild;
234}
235
236
237void VspInterior::SetupChildLinks(VspNode *b, VspNode *f)
238{
239    mBack = b;
240    mFront = f;
241}
242
243
244AxisAlignedBox3 VspInterior::GetBoundingBox() const
245{
246        return mBoundingBox;
247}
248
249
250void VspInterior::SetBoundingBox(const AxisAlignedBox3 &box)
251{
252        mBoundingBox = box;
253}
254
255
256int VspInterior::Type() const
257{
258        return Interior;
259}
260
261
262
263/****************************************************************/
264/*                  class VspLeaf implementation                */
265/****************************************************************/
266
267
268VspLeaf::VspLeaf(): mViewCell(NULL), mPvs(NULL)
269{
270}
271
272
273VspLeaf::~VspLeaf()
274{
275        DEL_PTR(mPvs);
276        CLEAR_CONTAINER(mVssRays);
277}
278
279
280int VspLeaf::Type() const
281{
282        return Leaf;
283}
284
285
286VspLeaf::VspLeaf(ViewCellLeaf *viewCell):
287mViewCell(viewCell)
288{
289}
290
291
292VspLeaf::VspLeaf(VspInterior *parent):
293VspNode(parent), mViewCell(NULL), mPvs(NULL)
294{}
295
296
297
298VspLeaf::VspLeaf(VspInterior *parent, ViewCellLeaf *viewCell):
299VspNode(parent), mViewCell(viewCell), mPvs(NULL)
300{
301}
302
303ViewCellLeaf *VspLeaf::GetViewCell() const
304{
305        return mViewCell;
306}
307
308void VspLeaf::SetViewCell(ViewCellLeaf *viewCell)
309{
310        mViewCell = viewCell;
311}
312
313
314bool VspLeaf::IsLeaf() const
315{
316        return true;
317}
318
319
320/******************************************************************************/
321/*                       class VspTree implementation                      */
322/******************************************************************************/
323
324
325VspTree::VspTree():
326mRoot(NULL),
327mOutOfBoundsCell(NULL),
328mStoreRays(false),
329mTimeStamp(1)
330{
331        bool randomize = false;
332        Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
333        if (randomize)
334                Randomize(); // initialise random generator for heuristics
335
336        //-- termination criteria for autopartition
337        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxDepth", mTermMaxDepth);
338        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minPvs", mTermMinPvs);
339        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minRays", mTermMinRays);
340        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minProbability", mTermMinProbability);
341        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxRayContribution", mTermMaxRayContribution);
342       
343        Environment::GetSingleton()->GetIntValue("VspTree.Termination.missTolerance", mTermMissTolerance);
344        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxViewCells", mMaxViewCells);
345
346        //-- max cost ratio for early tree termination
347        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxCostRatio", mTermMaxCostRatio);
348
349        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
350        Environment::GetSingleton()->GetIntValue("VspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
351
352        //-- factors for bsp tree split plane heuristics
353        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.ct_div_ci", mCtDivCi);
354
355        //-- partition criteria
356        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.epsilon", mEpsilon);
357
358        // if only the driving axis is used for axis aligned split
359        Environment::GetSingleton()->GetBoolValue("VspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
360       
361        //Environment::GetSingleton()->GetFloatValue("VspTree.maxTotalMemory", mMaxTotalMemory);
362        Environment::GetSingleton()->GetFloatValue("VspTree.maxStaticMemory", mMaxMemory);
363
364        Environment::GetSingleton()->GetBoolValue("VspTree.useCostHeuristics", mUseCostHeuristics);
365        Environment::GetSingleton()->GetBoolValue("VspTree.simulateOctree", mCirculatingAxis);
366       
367        Environment::GetSingleton()->GetIntValue("VspTree.pvsCountMethod", mPvsCountMethod);
368
369        char subdivisionStatsLog[100];
370        Environment::GetSingleton()->GetStringValue("VspTree.subdivisionStats", subdivisionStatsLog);
371        mSubdivisionStats.open(subdivisionStatsLog);
372
373        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.minBand", mMinBand);
374        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.maxBand", mMaxBand);
375       
376
377        //-- debug output
378
379        Debug << "******* VSP options ******** " << endl;
380    Debug << "max depth: " << mTermMaxDepth << endl;
381        Debug << "min PVS: " << mTermMinPvs << endl;
382        Debug << "min probabiliy: " << mTermMinProbability << endl;
383        Debug << "min rays: " << mTermMinRays << endl;
384        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
385        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
386        Debug << "miss tolerance: " << mTermMissTolerance << endl;
387        Debug << "max view cells: " << mMaxViewCells << endl;
388        Debug << "randomize: " << randomize << endl;
389
390        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
391        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
392        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
393        Debug << "max memory: " << mMaxMemory << endl;
394        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
395        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
396       
397        Debug << "circulating axis: " << mCirculatingAxis << endl;
398        Debug << "minband: " << mMinBand << endl;
399        Debug << "maxband: " << mMaxBand << endl;
400        Debug << "pvs count method: " << mPvsCountMethod << endl;
401       
402
403        mSplitCandidates = new vector<SortableEntry>;
404
405        Debug << endl;
406}
407
408
409VspViewCell *VspTree::GetOutOfBoundsCell()
410{
411        return mOutOfBoundsCell;
412}
413
414
415VspViewCell *VspTree::GetOrCreateOutOfBoundsCell()
416{
417        if (!mOutOfBoundsCell)
418        {
419                mOutOfBoundsCell = new VspViewCell();
420                mOutOfBoundsCell->SetId(-1);
421                mOutOfBoundsCell->SetValid(false);
422        }
423
424        return mOutOfBoundsCell;
425}
426
427
428const VspTreeStatistics &VspTree::GetStatistics() const
429{
430        return mVspStats;
431}
432
433
434VspTree::~VspTree()
435{
436        DEL_PTR(mRoot);
437        DEL_PTR(mSplitCandidates);
438}
439
440
441void VspTree::PrepareConstruction(const VssRayContainer &sampleRays,
442                                                                  AxisAlignedBox3 *forcedBoundingBox)
443{
444        mVspStats.nodes = 1;
445       
446        if (forcedBoundingBox)
447        {
448                mBoundingBox = *forcedBoundingBox;
449        }
450        else // compute vsp tree bounding box
451        {
452                mBoundingBox.Initialize();
453
454                VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
455
456                //-- compute bounding box
457        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
458                {
459                        VssRay *ray = *rit;
460
461                        // compute bounding box of view space
462                        mBoundingBox.Include(ray->GetTermination());
463                        mBoundingBox.Include(ray->GetOrigin());
464                }
465
466                mTermMinProbability *= mBoundingBox.GetVolume();
467                mGlobalCostMisses = 0;
468        }
469}
470
471
472void VspTree::AddSubdivisionStats(const int viewCells,
473                                                                  const float renderCostDecr,
474                                                                  const float splitCandidateCost,
475                                                                  const float totalRenderCost,
476                                                                  const float avgRenderCost)
477{
478        mSubdivisionStats
479                        << "#ViewCells\n" << viewCells << endl
480                        << "#RenderCostDecrease\n" << renderCostDecr << endl
481                        << "#SplitCandidateCost\n" << splitCandidateCost << endl
482                        << "#TotalRenderCost\n" << totalRenderCost << endl
483                        << "#AvgRenderCost\n" << avgRenderCost << endl;
484}
485
486
487// TODO: return memory usage in MB
488float VspTree::GetMemUsage() const
489{
490        return (float)
491                 (sizeof(VspTree) +
492                  mVspStats.Leaves() * sizeof(VspLeaf) +
493                  mCreatedViewCells * sizeof(VspViewCell) +
494                  mVspStats.pvs * sizeof(ObjectPvsData) +
495                  mVspStats.Interior() * sizeof(VspInterior) +
496                  mVspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
497}
498
499
500bool VspTree::LocalTerminationCriteriaMet(const VspTraversalData &data) const
501{
502        return
503#if TODO
504                (((int)data.mRays->size() <= mTermMinRays) ||
505                 (data.mPvs <= mTermMinPvs)   ||
506                 (data.mProbability <= mTermMinProbability) ||
507                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
508                 (data.mDepth >= mTermMaxDepth));
509#else
510                false;
511#endif
512}
513
514
515bool VspTree::GlobalTerminationCriteriaMet(const VspTraversalData &data) const
516{
517        return
518#if TODO
519                (mOutOfMemory ||
520                (mVspStats.Leaves() >= mMaxViewCells) ||
521                (mGlobalCostMisses >= mTermGlobalCostMissTolerance));
522#else
523                (mVspStats.Leaves() >= mMaxViewCells) ;         
524#endif
525}
526
527
528VspNode *VspTree::Subdivide(SplitQueue &tQueue,
529                                                        VspSplitCandidate &splitCandidate,
530                                                        const bool globalCriteriaMet)
531{
532        VspTraversalData &tData = splitCandidate.mParentData;
533
534        VspNode *newNode = tData.mNode;
535
536        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
537        {       
538                VspTraversalData tFrontData;
539                VspTraversalData tBackData;
540
541                //-- continue subdivision
542               
543                // create new interior node and two leaf node
544                const AxisAlignedPlane splitPlane = splitCandidate.mSplitPlane;
545                newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData);
546       
547                const int maxCostMisses = splitCandidate.mMaxCostMisses;
548
549
550                // how often was max cost ratio missed in this branch?
551                tFrontData.mMaxCostMisses = maxCostMisses;
552                tBackData.mMaxCostMisses = maxCostMisses;
553                       
554                //-- statistics
555                if (1)
556                {
557                        const float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
558                        const float cBack = (float)tBackData.mPvs * tBackData.mProbability;
559                        const float cData = (float)tData.mPvs * tData.mProbability;
560       
561                        const float costDecr =
562                                (cFront + cBack - cData) / mBoundingBox.GetVolume();
563
564                        mTotalCost += costDecr;
565                        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
566
567                        AddSubdivisionStats(mVspStats.Leaves(),
568                                                                -costDecr,
569                                                                splitCandidate.GetPriority(),
570                                                                mTotalCost,
571                                                                (float)mTotalPvsSize / (float)mVspStats.Leaves());
572                }
573
574       
575                //-- push the new split candidates on the queue
576                VspSplitCandidate *frontCandidate = new VspSplitCandidate();
577                VspSplitCandidate *backCandidate = new VspSplitCandidate();
578
579                EvalSplitCandidate(tFrontData, *frontCandidate);
580                EvalSplitCandidate(tBackData, *backCandidate);
581
582                tQueue.push(frontCandidate);
583                tQueue.push(backCandidate);
584
585                // delete old leaf node
586                DEL_PTR(tData.mNode);
587        }
588
589
590        //-- terminate traversal and create new view cell
591        if (newNode->IsLeaf())
592        {
593                VspLeaf *leaf = dynamic_cast<VspLeaf *>(newNode);
594       
595                VspViewCell *viewCell = new VspViewCell();
596        leaf->SetViewCell(viewCell);
597               
598                //-- update pvs
599                int conSamp = 0;
600                float sampCon = 0.0f;
601                AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
602
603                // update scalar pvs size value
604                mViewCellsManager->SetScalarPvsSize(viewCell, viewCell->GetPvs().GetSize());
605
606                mVspStats.contributingSamples += conSamp;
607                mVspStats.sampleContributions +=(int) sampCon;
608
609                //-- store additional info
610                if (mStoreRays)
611                {
612                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
613                        for (it = tData.mRays->begin(); it != it_end; ++ it)
614                        {
615                                (*it).mRay->Ref();                     
616                                leaf->mVssRays.push_back((*it).mRay);
617                        }
618                }
619               
620                viewCell->mLeaf = leaf;
621
622                viewCell->SetVolume(tData.mProbability);
623        leaf->mProbability = tData.mProbability;
624
625                // finally evaluate stats until this leaf
626                EvaluateLeafStats(tData);
627        }
628
629        //-- cleanup
630        tData.Clear();
631       
632        return newNode;
633}
634
635
636void VspTree::EvalSplitCandidate(VspTraversalData &tData,
637                                                                 VspSplitCandidate &splitCandidate)
638{
639        float frontProb;
640        float backProb;
641       
642        VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
643       
644        // compute locally best split plane
645        const bool success =
646                SelectSplitPlane(tData, splitCandidate.mSplitPlane,     frontProb, backProb);
647
648        // compute global decrease in render cost
649        splitCandidate.mPriority = EvalRenderCostDecrease(splitCandidate.mSplitPlane, tData);
650        splitCandidate.mParentData = tData;
651        splitCandidate.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
652
653        //Debug << "p: " << tData.mNode << " depth: " << tData.mDepth << endl;
654}
655
656
657VspInterior *VspTree::SubdivideNode(const AxisAlignedPlane &splitPlane,
658                                                                        VspTraversalData &tData,
659                                                                        VspTraversalData &frontData,
660                                                                        VspTraversalData &backData)
661{
662        VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
663       
664        //-- the front and back traversal data is filled with the new values
665        frontData.mDepth = tData.mDepth + 1;
666        frontData.mRays = new RayInfoContainer();
667       
668        backData.mDepth = tData.mDepth + 1;
669        backData.mRays = new RayInfoContainer();
670       
671        //-- subdivide rays
672        SplitRays(splitPlane,
673                          *tData.mRays,
674                          *frontData.mRays,
675                          *backData.mRays);
676
677        //Debug << "f: " << frontData.mRays->size() << " b: " << backData.mRays->size() << "d: " << tData.mRays->size() << endl;
678        //-- compute pvs
679        frontData.mPvs = ComputePvsSize(*frontData.mRays);
680        backData.mPvs = ComputePvsSize(*backData.mRays);
681
682        // split front and back node geometry and compute area
683        tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition,
684                                                         frontData.mBoundingBox, backData.mBoundingBox);
685
686
687        frontData.mProbability = frontData.mBoundingBox.GetVolume();
688        backData.mProbability = tData.mProbability - frontData.mProbability;
689
690       
691    ///////////////////////////////////////////
692        // subdivide further
693
694        // store maximal and minimal depth
695        if (tData.mDepth > mVspStats.maxDepth)
696        {
697                Debug << "max depth increases to " << tData.mDepth << " at " << mVspStats.Leaves() << " leaves" << endl;
698                mVspStats.maxDepth = tData.mDepth;
699        }
700
701        mVspStats.nodes += 2;
702
703   
704        VspInterior *interior = new VspInterior(splitPlane);
705
706#ifdef _DEBUG
707        Debug << interior << endl;
708#endif
709
710
711        //-- create front and back leaf
712
713        VspInterior *parent = leaf->GetParent();
714
715        // replace a link from node's parent
716        if (parent)
717        {
718                parent->ReplaceChildLink(leaf, interior);
719                interior->SetParent(parent);
720        }
721        else // new root
722        {
723                mRoot = interior;
724        }
725
726        // and setup child links
727        interior->SetupChildLinks(new VspLeaf(interior), new VspLeaf(interior));
728        // add bounding box
729        interior->SetBoundingBox(tData.mBoundingBox);
730
731        frontData.mNode = interior->GetFront();
732        backData.mNode = interior->GetBack();
733
734        interior->mTimeStamp = mTimeStamp ++;
735       
736        return interior;
737}
738
739
740void VspTree::AddToPvs(VspLeaf *leaf,
741                                                  const RayInfoContainer &rays,
742                                                  float &sampleContributions,
743                                                  int &contributingSamples)
744{
745        sampleContributions = 0;
746        contributingSamples = 0;
747 
748        RayInfoContainer::const_iterator it, it_end = rays.end();
749 
750        ViewCellLeaf *vc = leaf->GetViewCell();
751 
752        // add contributions from samples to the PVS
753        for (it = rays.begin(); it != it_end; ++ it)
754        {
755                float sc = 0.0f;
756                VssRay *ray = (*it).mRay;
757
758                bool madeContrib = false;
759                float contribution;
760
761                if (ray->mTerminationObject)
762                {
763                        if (vc->AddPvsSample(ray->mTerminationObject, ray->mPdf, contribution))
764                                madeContrib = true;
765                        sc += contribution;
766                }
767         
768                if (ray->mOriginObject)
769                {
770                        if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
771                                madeContrib = true;
772                        sc += contribution;
773                }
774         
775                sampleContributions += sc;
776
777                if (madeContrib)
778                        ++ contributingSamples;
779               
780                // store rays for visualization
781                if (0) leaf->mVssRays.push_back(new VssRay(*ray));
782        }
783}
784
785
786void VspTree::SortSplitCandidates(const RayInfoContainer &rays,
787                                                                  const int axis,
788                                                                  float minBand,
789                                                                  float maxBand)
790{
791        mSplitCandidates->clear();
792
793        int requestedSize = 2 * (int)(rays.size());
794
795        // creates a sorted split candidates array
796        if (mSplitCandidates->capacity() > 500000 &&
797                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
798        {
799        delete mSplitCandidates;
800                mSplitCandidates = new vector<SortableEntry>;
801        }
802
803        mSplitCandidates->reserve(requestedSize);
804
805        float pos;
806
807        // insert all queries
808        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
809        {
810                const bool positive = (*ri).mRay->HasPosDir(axis);
811                               
812                pos = (*ri).ExtrapOrigin(axis);
813               
814                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
815                                                                        pos, (*ri).mRay));
816
817                pos = (*ri).ExtrapTermination(axis);
818
819                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
820                                                                        pos, (*ri).mRay));
821        }
822
823        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
824}
825
826
827int VspTree::GetPvsContribution(Intersectable *object) const
828{
829    int pvsContri = 0;
830
831        KdPvsMap::const_iterator kit, kit_end = object->mKdPvs.mEntries.end();
832
833        Intersectable::NewMail();
834
835        // Search kd leaves this object is attached to
836        for (kit = object->mKdPvs.mEntries.begin(); kit != kit_end; ++ kit)
837        {
838                KdNode *l = (*kit).first;
839
840                // new object found during sweep
841                // => increase pvs contribution of this kd node
842                if (!l->Mailed())
843                {
844                        l->Mail();
845                        ++ pvsContri;
846                }
847        }
848
849        return pvsContri;
850}
851
852
853int VspTree::PrepareHeuristics(Intersectable *object)
854{
855        set<KdLeaf *>::const_iterator kit, kit_end = object->mKdLeaves.end();
856       
857        int pvsSize = 0;
858
859        for (kit = object->mKdLeaves.begin(); kit != kit_end; ++ kit)
860        {
861                KdLeaf *node = dynamic_cast<KdLeaf *>(*kit);
862
863                if (!node->Mailed())
864                {
865                        node->Mail();
866                        node->mCounter = 1;
867
868                        //Debug << "here5 "<<(int)node->mObjects.size() <<" "<< node->mMultipleObjects.size()<< " " <<(int)node->mObjects.size() - node->mMultipleObjects.size()<<endl;
869                        // add objects without the objects which are in several kd leaves
870                        pvsSize += (int)(node->mObjects.size() - node->mMultipleObjects.size());
871                }
872                else
873                {
874                        ++ node->mCounter;
875                }
876
877                //-- the objects belonging to several leaves must be handled seperately
878                ObjectContainer::const_iterator oit, oit_end = node->mMultipleObjects.end();
879
880                for (oit = node->mMultipleObjects.begin(); oit != oit_end; ++ oit)
881                {
882                        Intersectable *object = *oit;
883                                               
884                        if (!object->Mailed())
885                        {
886                                //Debug << "here233: " << object->mKdLeaves.size() << endl;
887                                object->Mail();
888                                object->mCounter = 1;
889
890                                ++ pvsSize;
891                        }
892                        else
893                        {
894                                ++ object->mCounter;
895                        }
896                }
897        }
898
899        return pvsSize;
900}
901
902
903int VspTree::PrepareHeuristics(const RayInfoContainer &rays)
904{       
905        Intersectable::NewMail();
906        KdNode::NewMail();
907
908        int pvsSize = 0;
909
910        RayInfoContainer::const_iterator ri, ri_end = rays.end();
911
912        //-- set all kd nodes as belonging to the front pvs
913
914        for (ri = rays.begin(); ri != ri_end; ++ ri)
915        {
916                Intersectable *oObject = (*ri).mRay->mOriginObject;
917
918                if (oObject)
919                {
920                        if (mPvsCountMethod == PER_OBJECT)
921                        {
922                                if (!oObject->Mailed())
923                                {
924                                        oObject->Mail();
925                                        oObject->mCounter = 1;
926                                        ++ pvsSize;
927                                }
928                                else
929                                {
930                                        ++ oObject->mCounter;
931                                }
932                        }
933                        else
934                        {
935                                pvsSize += PrepareHeuristics(oObject); 
936                        }       
937                }
938
939                Intersectable *tObject = (*ri).mRay->mTerminationObject;
940
941                if (tObject)
942                {
943                        if (mPvsCountMethod == PER_OBJECT)
944                        {
945                                if (!tObject->Mailed())
946                                {
947                                        tObject->Mail();
948                                        tObject->mCounter = 1;
949                                        ++ pvsSize;
950                                }
951                                else
952                                {
953                                        ++ tObject->mCounter;
954                                }
955                        }
956                        else
957                        {
958                                pvsSize += PrepareHeuristics(tObject); 
959                        }       
960                }
961        }
962
963        return pvsSize;
964}
965
966
967int VspTree::GetPvsIncr(Intersectable *object, const KdPvsMap &activeNodes)
968{
969        // TODO:
970        // use kd info table to apply set theory in order to
971        // sort out dublicate pvs entries due to objects which
972        // belong to more than one kd leaves.
973#if TODO
974        KdPvsMap::const_iterator kit, kit_end = obj->mKdPvs.mEntries.end();
975
976        // Search kd leaves this object is attached to
977        for (kit = obj->mKdPvs.mEntries.begin(); kit != kit_end; ++ kit)
978        {
979                KdNode *l = (*kit).first;
980
981                // new object found during sweep
982                // => increase pvs contribution of this kd node
983                if (!l->Mailed())
984                {
985                        l->Mail();
986                        ++ pvsContr;
987                }
988        }
989#endif
990        return 0;
991}
992
993
994void VspTree::RemoveContriFromPvs(KdLeaf *leaf, int &pvs) const
995{
996        // leaf falls out of right pvs
997        if (-- leaf->mCounter == 0)
998        {
999                pvs -= ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
1000        }
1001
1002        //-- handle objects which are in several kd leaves separately
1003        ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1004
1005        for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1006        {
1007                Intersectable *object = *oit;
1008
1009                if (-- object->mCounter == 0)
1010                {
1011                        -- pvs;
1012                }
1013        }
1014}
1015
1016
1017void VspTree::AddContriToPvs(KdLeaf *leaf, int &pvs) const
1018{
1019        if (!leaf->Mailed())
1020        {
1021                leaf->Mail();
1022
1023                // add objects without those which are part of several kd leaves
1024                pvs += ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
1025
1026                //-- handle objects of several kd leaves separately
1027                ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1028
1029                for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1030                {
1031                        Intersectable *object = *oit;
1032
1033                        // object not previously in left pvs
1034                        if (!object->Mailed())
1035                        {
1036                                object->Mail();
1037                                ++ pvs;
1038                        }
1039                }
1040        }                                       
1041}
1042
1043
1044void VspTree::EvalPvsIncr(const SortableEntry &ci,
1045                                                  int &pvsLeft,
1046                                                  int &pvsRight) const
1047{
1048        VssRay *ray = ci.ray;
1049
1050        Intersectable *oObject = ray->mOriginObject;
1051        Intersectable *tObject = ray->mTerminationObject;
1052
1053        if (oObject)
1054        {       
1055                if (mPvsCountMethod == PER_OBJECT)
1056                {
1057                        if (ci.type == SortableEntry::ERayMin)
1058                        {
1059                                if (!oObject->Mailed())
1060                                {
1061                                        oObject->Mail();
1062                                        ++ pvsLeft;
1063                                }
1064                        }
1065                        else if (ci.type == SortableEntry::ERayMax)
1066                        {
1067                                if (-- oObject->mCounter == 0)
1068                                        -- pvsRight;
1069                        }
1070                }
1071                else // per kd node
1072                {
1073                        set<KdLeaf *>::const_iterator kit, kit_end = oObject->mKdLeaves.end();
1074
1075                        for (kit = oObject->mKdLeaves.begin(); kit != kit_end; ++ kit)
1076                        {
1077                                KdLeaf *node = dynamic_cast<KdLeaf *>(*kit);
1078
1079                                // add contributions of the kd nodes
1080                                if (ci.type == SortableEntry::ERayMin)
1081                                {
1082                                        AddContriToPvs(node, pvsLeft);
1083                                }
1084                                else if (ci.type == SortableEntry::ERayMax)
1085                                {
1086                                        RemoveContriFromPvs(node, pvsRight);
1087                                }
1088                        }
1089                }
1090                       
1091
1092        }
1093       
1094        if (tObject)
1095        {       
1096                if (mPvsCountMethod == PER_OBJECT)
1097                {
1098                        if (ci.type == SortableEntry::ERayMin)
1099                        {
1100                                if (!tObject->Mailed())
1101                                {
1102                                        tObject->Mail();
1103                                        ++ pvsLeft;
1104                                }
1105                        }
1106                        else if (ci.type == SortableEntry::ERayMax)
1107                        {
1108                                if (-- tObject->mCounter == 0)
1109                                        -- pvsRight;
1110                        }
1111                }
1112                else // per kd node
1113                {
1114                        set<KdLeaf *>::const_iterator kit, kit_end = tObject->mKdLeaves.end();
1115                       
1116                        for (kit = tObject->mKdLeaves.begin(); kit != kit_end; ++ kit)
1117                        {
1118                                KdLeaf *node = dynamic_cast<KdLeaf *>(*kit);
1119
1120                                if (ci.type == SortableEntry::ERayMin)
1121                                {
1122                                        AddContriToPvs(node, pvsLeft);
1123                                }
1124                                else if (ci.type == SortableEntry::ERayMax)
1125                                {
1126                                        RemoveContriFromPvs(node, pvsRight);
1127                                }
1128                        }
1129                }
1130        }
1131}
1132
1133
1134float VspTree::EvalLocalCostHeuristics(const RayInfoContainer &rays,
1135                                                                           const AxisAlignedBox3 &box,
1136                                                                           int pvsSize,
1137                                                                           const int axis,
1138                                                                           float &position)
1139{
1140        const float minBox = box.Min(axis);
1141        const float maxBox = box.Max(axis);
1142
1143        const float sizeBox = maxBox - minBox;
1144
1145        const float minBand = minBox + mMinBand * sizeBox;
1146        const float maxBand = minBox + mMaxBand * sizeBox;
1147
1148        SortSplitCandidates(rays, axis, minBand, maxBand);
1149
1150        // prepare the sweep
1151        // (note: returns pvs size, so there would be no need
1152        // to give pvs size as argument)
1153        pvsSize = PrepareHeuristics(rays);
1154
1155        Debug << "here45 pvs: " << pvsSize << endl;
1156
1157        // go through the lists, count the number of objects left and right
1158        // and evaluate the following cost funcion:
1159        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
1160
1161        int pvsl = 0;
1162        int pvsr = pvsSize;
1163
1164        int pvsBack = pvsl;
1165        int pvsFront = pvsr;
1166
1167        float sum = (float)pvsSize * sizeBox;
1168        float minSum = 1e20f;
1169
1170       
1171        // if no good split can be found, take mid split
1172        position = minBox + 0.5f * sizeBox;
1173       
1174        // the relative cost ratio
1175        float ratio = 99999999.0f;
1176        bool splitPlaneFound = false;
1177
1178        Intersectable::NewMail();
1179        KdLeaf::NewMail();
1180
1181
1182        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
1183
1184Debug << "****************" << endl;
1185        //-- traverse through visibility events
1186
1187        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
1188        {
1189                EvalPvsIncr(*ci, pvsl, pvsr);
1190
1191                // Note: sufficient to compare size of bounding boxes of front and back side?
1192                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
1193                {
1194                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
1195
1196                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << "\t cost= " << sum << endl;
1197
1198                        if (sum < minSum)
1199                        {
1200                                splitPlaneFound = true;
1201
1202                                minSum = sum;
1203                                position = (*ci).value;
1204                               
1205                                pvsBack = pvsl;
1206                                pvsFront = pvsr;
1207                        }
1208                }
1209        }
1210       
1211       
1212        // -- compute cost
1213
1214        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1215        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1216
1217        const float pOverall = sizeBox;
1218        const float pBack = position - minBox;
1219        const float pFront = maxBox - position;
1220
1221        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
1222    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1223        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1224       
1225        const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
1226        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1227
1228        if (splitPlaneFound)
1229        {
1230                ratio = newRenderCost / oldRenderCost;
1231        }
1232       
1233        //if (axis != 1)
1234        Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1235              <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
1236
1237        return ratio;
1238}
1239
1240
1241float VspTree::SelectSplitPlane(const VspTraversalData &tData,
1242                                                                AxisAlignedPlane &plane,
1243                                                                float &pFront,
1244                                                                float &pBack)
1245{
1246        float nPosition[3];
1247        float nCostRatio[3];
1248        float nProbFront[3];
1249        float nProbBack[3];
1250
1251        // create bounding box of node geometry
1252        AxisAlignedBox3 box = tData.mBoundingBox;
1253               
1254        int sAxis = 0;
1255        int bestAxis = -1;
1256
1257        // if we use some kind of specialised fixed axis
1258    const bool useSpecialAxis =
1259                mOnlyDrivingAxis || mCirculatingAxis;
1260        //Debug << "data: " << tData.mBoundingBox << " pvs " << tData.mPvs << endl;
1261        if (mCirculatingAxis)
1262        {
1263                int parentAxis = 0;
1264                VspNode *parent = tData.mNode->GetParent();
1265
1266                if (parent)
1267                        parentAxis = dynamic_cast<VspInterior *>(parent)->GetAxis();
1268
1269                sAxis = (parentAxis + 1) % 3;
1270        }
1271        else if (mOnlyDrivingAxis)
1272        {
1273                sAxis = box.Size().DrivingAxis();
1274        }
1275        //sAxis = 2;
1276        for (int axis = 0; axis < 3; ++ axis)
1277        {
1278                if (!useSpecialAxis || (axis == sAxis))
1279                {
1280                        //-- place split plane using heuristics
1281
1282                        if (mUseCostHeuristics)
1283                        {
1284                                nCostRatio[axis] =
1285                                        EvalLocalCostHeuristics(*tData.mRays,
1286                                                                                        box,
1287                                                                                        tData.mPvs,
1288                                                                                        axis,
1289                                                                                        nPosition[axis]);                       
1290                        }
1291                        else //-- split plane position is spatial median
1292                        {
1293                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1294
1295                                nCostRatio[axis] = EvalLocalSplitCost(tData,
1296                                                                                                          box,
1297                                                                                                          axis,
1298                                                                                                          nPosition[axis],
1299                                                                                                          nProbFront[axis],
1300                                                                                                          nProbBack[axis]);
1301                        }
1302                                               
1303                        if (bestAxis == -1)
1304                        {
1305                                bestAxis = axis;
1306                        }
1307                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1308                        {
1309                                bestAxis = axis;
1310                        }
1311                }
1312        }
1313
1314
1315        //-- assign values
1316       
1317        plane.mAxis = bestAxis;
1318        // split plane position
1319    plane.mPosition = nPosition[bestAxis];
1320
1321        pFront = nProbFront[bestAxis];
1322        pBack = nProbBack[bestAxis];
1323
1324        //Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
1325        return nCostRatio[bestAxis];
1326}
1327
1328
1329inline void VspTree::GenerateUniqueIdsForPvs()
1330{
1331        Intersectable::NewMail(); sBackId = Intersectable::sMailId;
1332        Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
1333        Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
1334}
1335
1336
1337float VspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
1338                                                                          const VspTraversalData &data) const
1339{
1340#if 0
1341        return (float)-data.mDepth;
1342#endif
1343        float pvsFront = 0;
1344        float pvsBack = 0;
1345        float totalPvs = 0;
1346
1347        // probability that view point lies in back / front node
1348        float pOverall = data.mProbability;
1349        float pFront = 0;
1350        float pBack = 0;
1351
1352
1353        // create unique ids for pvs heuristics
1354        GenerateUniqueIdsForPvs();
1355       
1356        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1357
1358        for (rit = data.mRays->begin(); rit != rit_end; ++ rit)
1359        {
1360                RayInfo rayInf = *rit;
1361
1362                float t;
1363                VssRay *ray = rayInf.mRay;
1364
1365                const int cf =
1366                        rayInf.ComputeRayIntersection(candidatePlane.mAxis,
1367                                                                                  candidatePlane.mPosition, t);
1368
1369                // find front and back pvs for origing and termination object
1370                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1371                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1372        }
1373
1374
1375        AxisAlignedBox3 frontBox;
1376        AxisAlignedBox3 backBox;
1377
1378        data.mBoundingBox.Split(candidatePlane.mAxis, candidatePlane.mPosition, frontBox, backBox);
1379
1380        pFront = frontBox.GetVolume();
1381        pBack = pOverall - pFront;
1382               
1383
1384        //-- pvs rendering heuristics
1385        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1386        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1387
1388        //-- only render cost heuristics or combined with standard deviation
1389        const float penaltyOld = EvalPvsPenalty((int)totalPvs, lowerPvsLimit, upperPvsLimit);
1390    const float penaltyFront = EvalPvsPenalty((int)pvsFront, lowerPvsLimit, upperPvsLimit);
1391        const float penaltyBack = EvalPvsPenalty((int)pvsBack, lowerPvsLimit, upperPvsLimit);
1392                       
1393        const float oldRenderCost = pOverall * penaltyOld;
1394        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1395
1396        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
1397        const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBoundingBox.GetVolume();
1398       
1399        // take render cost of node into account
1400        // otherwise danger of being stuck in a local minimum!!
1401        const float factor = 0.99f;
1402
1403        const float normalizedOldRenderCost = oldRenderCost / mBoundingBox.GetVolume();
1404        return factor * renderCostDecrease + (1.0f - factor) * normalizedOldRenderCost;
1405}
1406
1407
1408float VspTree::EvalLocalSplitCost(const VspTraversalData &data,
1409                                                                  const AxisAlignedBox3 &box,
1410                                                                  const int axis,
1411                                                                  const float &position,
1412                                                                  float &pFront,
1413                                                                  float &pBack) const
1414{
1415        float pvsTotal = 0;
1416        float pvsFront = 0;
1417        float pvsBack = 0;
1418       
1419        // create unique ids for pvs heuristics
1420        GenerateUniqueIdsForPvs();
1421
1422        const int pvsSize = data.mPvs;
1423
1424        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1425
1426        // this is the main ray classification loop!
1427        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
1428        {
1429                // determine the side of this ray with respect to the plane
1430                float t;
1431                const int side = (*rit).ComputeRayIntersection(axis, position, t);
1432       
1433                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
1434                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
1435        }
1436
1437        //-- pvs heuristics
1438        float pOverall;
1439
1440        //-- compute heurstics
1441       
1442        pOverall = data.mProbability;
1443        // we take simplified computation for mid split
1444        pBack = pFront = pOverall * 0.5f;
1445       
1446       
1447#ifdef _DEBUG
1448        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
1449        Debug << pFront << " " << pBack << " " << pOverall << endl;
1450#endif
1451
1452       
1453        const float newCost = pvsBack * pBack + pvsFront * pFront;
1454        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
1455
1456        return  (mCtDivCi + newCost) / oldCost;
1457}
1458
1459
1460void VspTree::AddObjToPvs(Intersectable *obj,
1461                                                  const int cf,
1462                                                  float &frontPvs,
1463                                                  float &backPvs,
1464                                                  float &totalPvs) const
1465{
1466        if (!obj)
1467                return;
1468
1469        //const float renderCost = mViewCellsManager->EvalRenderCost(obj);
1470        const int renderCost = 1;
1471
1472        // new object
1473        if ((obj->mMailbox != sFrontId) &&
1474                (obj->mMailbox != sBackId) &&
1475                (obj->mMailbox != sFrontAndBackId))
1476        {
1477                totalPvs += renderCost;
1478        }
1479
1480        // TODO: does this really belong to no pvs?
1481        //if (cf == Ray::COINCIDENT) return;
1482
1483        // object belongs to both PVS
1484        if (cf >= 0)
1485        {
1486                if ((obj->mMailbox != sFrontId) &&
1487                        (obj->mMailbox != sFrontAndBackId))
1488                {
1489                        frontPvs += renderCost;
1490               
1491                        if (obj->mMailbox == sBackId)
1492                                obj->mMailbox = sFrontAndBackId;
1493                        else
1494                                obj->mMailbox = sFrontId;
1495                }
1496        }
1497
1498        if (cf <= 0)
1499        {
1500                if ((obj->mMailbox != sBackId) &&
1501                        (obj->mMailbox != sFrontAndBackId))
1502                {
1503                        backPvs += renderCost;
1504               
1505                        if (obj->mMailbox == sFrontId)
1506                                obj->mMailbox = sFrontAndBackId;
1507                        else
1508                                obj->mMailbox = sBackId;
1509                }
1510        }
1511}
1512
1513
1514void VspTree::CollectLeaves(vector<VspLeaf *> &leaves,
1515                                                           const bool onlyUnmailed,
1516                                                           const int maxPvsSize) const
1517{
1518        stack<VspNode *> nodeStack;
1519        nodeStack.push(mRoot);
1520
1521        while (!nodeStack.empty())
1522        {
1523                VspNode *node = nodeStack.top();
1524                nodeStack.pop();
1525               
1526                if (node->IsLeaf())
1527                {
1528                        // test if this leaf is in valid view space
1529                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1530                        if (leaf->TreeValid() &&
1531                                (!onlyUnmailed || !leaf->Mailed()) &&
1532                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
1533                        {
1534                                leaves.push_back(leaf);
1535                        }
1536                }
1537                else
1538                {
1539                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1540
1541                        nodeStack.push(interior->GetBack());
1542                        nodeStack.push(interior->GetFront());
1543                }
1544        }
1545}
1546
1547
1548AxisAlignedBox3 VspTree::GetBoundingBox() const
1549{
1550        return mBoundingBox;
1551}
1552
1553
1554VspNode *VspTree::GetRoot() const
1555{
1556        return mRoot;
1557}
1558
1559
1560void VspTree::EvaluateLeafStats(const VspTraversalData &data)
1561{
1562        // the node became a leaf -> evaluate stats for leafs
1563        VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);
1564
1565
1566        if (data.mPvs > mVspStats.maxPvs)
1567        {
1568                mVspStats.maxPvs = data.mPvs;
1569        }
1570
1571        mVspStats.pvs += data.mPvs;
1572
1573        if (data.mDepth < mVspStats.minDepth)
1574        {
1575                mVspStats.minDepth = data.mDepth;
1576        }
1577       
1578        if (data.mDepth >= mTermMaxDepth)
1579        {
1580        ++ mVspStats.maxDepthNodes;
1581                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
1582        }
1583
1584        // accumulate rays to compute rays /  leaf
1585        mVspStats.accumRays += (int)data.mRays->size();
1586
1587        if (data.mPvs < mTermMinPvs)
1588                ++ mVspStats.minPvsNodes;
1589
1590        if ((int)data.mRays->size() < mTermMinRays)
1591                ++ mVspStats.minRaysNodes;
1592
1593        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1594                ++ mVspStats.maxRayContribNodes;
1595
1596        if (data.mProbability <= mTermMinProbability)
1597                ++ mVspStats.minProbabilityNodes;
1598       
1599        // accumulate depth to compute average depth
1600        mVspStats.accumDepth += data.mDepth;
1601
1602        ++ mCreatedViewCells;
1603
1604#ifdef _DEBUG
1605        Debug << "BSP stats: "
1606                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1607                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1608                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1609                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "), "
1610                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1611#endif
1612}
1613
1614
1615void VspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
1616{
1617        ViewCell::NewMail();
1618        CollectViewCells(mRoot, onlyValid, viewCells, true);
1619}
1620
1621
1622void VspTree::CollapseViewCells()
1623{
1624// TODO
1625#if HAS_TO_BE_REDONE
1626        stack<VspNode *> nodeStack;
1627
1628        if (!mRoot)
1629                return;
1630
1631        nodeStack.push(mRoot);
1632       
1633        while (!nodeStack.empty())
1634        {
1635                VspNode *node = nodeStack.top();
1636                nodeStack.pop();
1637               
1638                if (node->IsLeaf())
1639        {
1640                        BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1641
1642                        if (!viewCell->GetValid())
1643                        {
1644                                BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1645       
1646                                ViewCellContainer leaves;
1647                                mViewCellsTree->CollectLeaves(viewCell, leaves);
1648
1649                                ViewCellContainer::const_iterator it, it_end = leaves.end();
1650
1651                                for (it = leaves.begin(); it != it_end; ++ it)
1652                                {
1653                                        VspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
1654                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
1655                                        ++ mVspStats.invalidLeaves;
1656                                }
1657
1658                                // add to unbounded view cell
1659                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
1660                                DEL_PTR(viewCell);
1661                        }
1662                }
1663                else
1664                {
1665                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1666               
1667                        nodeStack.push(interior->GetFront());
1668                        nodeStack.push(interior->GetBack());
1669                }
1670        }
1671
1672        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1673#endif
1674}
1675
1676
1677void VspTree::CollectRays(VssRayContainer &rays)
1678{
1679        vector<VspLeaf *> leaves;
1680
1681        vector<VspLeaf *>::const_iterator lit, lit_end = leaves.end();
1682
1683        for (lit = leaves.begin(); lit != lit_end; ++ lit)
1684        {
1685                VspLeaf *leaf = *lit;
1686                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
1687
1688                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
1689                        rays.push_back(*rit);
1690        }
1691}
1692
1693
1694void VspTree::SetViewCellsManager(ViewCellsManager *vcm)
1695{
1696        mViewCellsManager = vcm;
1697}
1698
1699
1700void VspTree::ValidateTree()
1701{
1702        mVspStats.invalidLeaves = 0;
1703
1704        stack<VspNode *> nodeStack;
1705
1706        if (!mRoot)
1707                return;
1708
1709        nodeStack.push(mRoot);
1710
1711        while (!nodeStack.empty())
1712        {
1713                VspNode *node = nodeStack.top();
1714                nodeStack.pop();
1715               
1716                if (node->IsLeaf())
1717                {
1718                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1719
1720                        if (!leaf->GetViewCell()->GetValid())
1721                                ++ mVspStats.invalidLeaves;
1722
1723                        // validity flags don't match => repair
1724                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
1725                        {
1726                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
1727                                PropagateUpValidity(leaf);
1728                        }
1729                }
1730                else
1731                {
1732                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1733               
1734                        nodeStack.push(interior->GetFront());
1735                        nodeStack.push(interior->GetBack());
1736                }
1737        }
1738
1739        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1740}
1741
1742
1743
1744void VspTree::CollectViewCells(VspNode *root,
1745                                                                  bool onlyValid,
1746                                                                  ViewCellContainer &viewCells,
1747                                                                  bool onlyUnmailed) const
1748{
1749        stack<VspNode *> nodeStack;
1750
1751        if (!root)
1752                return;
1753
1754        nodeStack.push(root);
1755       
1756        while (!nodeStack.empty())
1757        {
1758                VspNode *node = nodeStack.top();
1759                nodeStack.pop();
1760               
1761                if (node->IsLeaf())
1762                {
1763                        if (!onlyValid || node->TreeValid())
1764                        {
1765                                ViewCellLeaf *leafVc = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1766
1767                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
1768                                               
1769                                if (!onlyUnmailed || !viewCell->Mailed())
1770                                {
1771                                        viewCell->Mail();
1772                                        viewCells.push_back(viewCell);
1773                                }
1774                        }
1775                }
1776                else
1777                {
1778                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1779               
1780                        nodeStack.push(interior->GetFront());
1781                        nodeStack.push(interior->GetBack());
1782                }
1783        }
1784}
1785
1786
1787int VspTree::FindNeighbors(VspLeaf *n,
1788                                                   vector<VspLeaf *> &neighbors,
1789                                                   const bool onlyUnmailed) const
1790{
1791        stack<VspNode *> nodeStack;
1792        nodeStack.push(mRoot);
1793
1794        const AxisAlignedBox3 box = GetBBox(n);
1795
1796        while (!nodeStack.empty())
1797        {
1798                VspNode *node = nodeStack.top();
1799                nodeStack.pop();
1800
1801                if (node->IsLeaf())
1802                {
1803                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1804
1805                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
1806                                neighbors.push_back(leaf);
1807                }
1808                else
1809                {
1810                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1811                       
1812                        if (interior->GetPosition() > box.Max(interior->GetAxis()))
1813                                nodeStack.push(interior->GetBack());
1814                        else
1815                        {
1816                                if (interior->GetPosition() < box.Min(interior->GetAxis()))
1817                                        nodeStack.push(interior->GetFront());
1818                                else
1819                                {
1820                                        // random decision
1821                                        nodeStack.push(interior->GetBack());
1822                                        nodeStack.push(interior->GetFront());
1823                                }
1824                        }
1825                }
1826        }
1827
1828        return (int)neighbors.size();
1829}
1830
1831
1832// Find random neighbor which was not mailed
1833VspLeaf *VspTree::GetRandomLeaf(const Plane3 &plane)
1834{
1835        stack<VspNode *> nodeStack;
1836        nodeStack.push(mRoot);
1837 
1838        int mask = rand();
1839 
1840        while (!nodeStack.empty())
1841        {
1842                VspNode *node = nodeStack.top();
1843               
1844                nodeStack.pop();
1845               
1846                if (node->IsLeaf())
1847                {
1848                        return dynamic_cast<VspLeaf *>(node);
1849                }
1850                else
1851                {
1852                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1853                        VspNode *next;
1854                       
1855                        if (GetBBox(interior->GetBack()).Side(plane) < 0)
1856                        {
1857                                next = interior->GetFront();
1858                        }
1859            else
1860                        {
1861                                if (GetBBox(interior->GetFront()).Side(plane) < 0)
1862                                {
1863                                        next = interior->GetBack();
1864                                }
1865                                else
1866                                {
1867                                        // random decision
1868                                        if (mask & 1)
1869                                                next = interior->GetBack();
1870                                        else
1871                                                next = interior->GetFront();
1872                                        mask = mask >> 1;
1873                                }
1874                        }
1875                       
1876                        nodeStack.push(next);
1877                }
1878        }
1879 
1880        return NULL;
1881}
1882
1883
1884VspLeaf *VspTree::GetRandomLeaf(const bool onlyUnmailed)
1885{
1886        stack<VspNode *> nodeStack;
1887
1888        nodeStack.push(mRoot);
1889
1890        int mask = rand();
1891
1892        while (!nodeStack.empty())
1893        {
1894                VspNode *node = nodeStack.top();
1895                nodeStack.pop();
1896
1897                if (node->IsLeaf())
1898                {
1899                        if ( (!onlyUnmailed || !node->Mailed()) )
1900                                return dynamic_cast<VspLeaf *>(node);
1901                }
1902                else
1903                {
1904                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1905
1906                        // random decision
1907                        if (mask & 1)
1908                                nodeStack.push(interior->GetBack());
1909                        else
1910                                nodeStack.push(interior->GetFront());
1911
1912                        mask = mask >> 1;
1913                }
1914        }
1915
1916        return NULL;
1917}
1918
1919
1920void VspTree::CollectPvs(const RayInfoContainer &rays,
1921                                                 ObjectContainer &objects) const
1922{
1923        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1924
1925        Intersectable::NewMail();
1926
1927        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1928        {
1929                VssRay *ray = (*rit).mRay;
1930
1931                Intersectable *object;
1932                object = ray->mOriginObject;
1933
1934        if (object)
1935                {
1936                        if (!object->Mailed())
1937                        {
1938                                object->Mail();
1939                                objects.push_back(object);
1940                        }
1941                }
1942
1943                object = ray->mTerminationObject;
1944
1945                if (object)
1946                {
1947                        if (!object->Mailed())
1948                        {
1949                                object->Mail();
1950                                objects.push_back(object);
1951                        }
1952                }
1953        }
1954}
1955
1956
1957int VspTree::ComputePvsSize(const RayInfoContainer &rays) const
1958{
1959        int pvsSize = 0;
1960
1961        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1962
1963        Intersectable::NewMail();
1964
1965        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1966        {
1967                VssRay *ray = (*rit).mRay;
1968
1969                if (ray->mOriginObject)
1970                {
1971                        if (!ray->mOriginObject->Mailed())
1972                        {
1973                                ray->mOriginObject->Mail();
1974                                ++ pvsSize;
1975                        }
1976                }
1977                if (ray->mTerminationObject)
1978                {
1979                        if (!ray->mTerminationObject->Mailed())
1980                        {
1981                                ray->mTerminationObject->Mail();
1982                                ++ pvsSize;
1983                        }
1984                }
1985        }
1986
1987        return pvsSize;
1988}
1989
1990
1991float VspTree::GetEpsilon() const
1992{
1993        return mEpsilon;
1994}
1995
1996
1997int VspTree::CastLineSegment(const Vector3 &origin,
1998                                                         const Vector3 &termination,
1999                             ViewCellContainer &viewcells)
2000{
2001        int hits = 0;
2002
2003        float mint = 0.0f, maxt = 1.0f;
2004        const Vector3 dir = termination - origin;
2005
2006        stack<LineTraversalData> tStack;
2007
2008        Intersectable::NewMail();
2009        ViewCell::NewMail();
2010
2011        Vector3 entp = origin;
2012        Vector3 extp = termination;
2013
2014        VspNode *node = mRoot;
2015        VspNode *farChild;
2016
2017        float position;
2018        int axis;
2019
2020        while (1)
2021        {
2022                if (!node->IsLeaf())
2023                {
2024                        VspInterior *in = dynamic_cast<VspInterior *>(node);
2025                        position = in->GetPosition();
2026                        axis = in->GetAxis();
2027
2028                        if (entp[axis] <= position)
2029                        {
2030                                if (extp[axis] <= position)
2031                                {
2032                                        node = in->GetBack();
2033                                        // cases N1,N2,N3,P5,Z2,Z3
2034                                        continue;
2035                                } else
2036                                {
2037                                        // case N4
2038                                        node = in->GetBack();
2039                                        farChild = in->GetFront();
2040                                }
2041                        }
2042                        else
2043                        {
2044                                if (position <= extp[axis])
2045                                {
2046                                        node = in->GetFront();
2047                                        // cases P1,P2,P3,N5,Z1
2048                                        continue;
2049                                }
2050                                else
2051                                {
2052                                        node = in->GetFront();
2053                                        farChild = in->GetBack();
2054                                        // case P4
2055                                }
2056                        }
2057
2058                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2059                        // case N4 or P4
2060                        const float tdist = (position - origin[axis]) / dir[axis];
2061                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2062
2063                        extp = origin + dir * tdist;
2064                        maxt = tdist;
2065                }
2066                else
2067                {
2068                        // compute intersection with all objects in this leaf
2069                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2070                        ViewCell *vc = leaf->GetViewCell();
2071
2072                        if (!vc->Mailed())
2073                        {
2074                                vc->Mail();
2075                                viewcells.push_back(vc);
2076                                ++ hits;
2077                        }
2078#if 0
2079                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2080#endif
2081                        // get the next node from the stack
2082                        if (tStack.empty())
2083                                break;
2084
2085                        entp = extp;
2086                        mint = maxt;
2087                       
2088                        LineTraversalData &s  = tStack.top();
2089                        node = s.mNode;
2090                        extp = s.mExitPoint;
2091                        maxt = s.mMaxT;
2092                        tStack.pop();
2093                }
2094        }
2095
2096        return hits;
2097}
2098
2099
2100int VspTree::CastRay(Ray &ray)
2101{
2102        int hits = 0;
2103
2104        stack<LineTraversalData> tStack;
2105        const Vector3 dir = ray.GetDir();
2106
2107        float maxt, mint;
2108
2109        if (!mBoundingBox.GetRaySegment(ray, mint, maxt))
2110                return 0;
2111
2112        Intersectable::NewMail();
2113        ViewCell::NewMail();
2114
2115        Vector3 entp = ray.Extrap(mint);
2116        Vector3 extp = ray.Extrap(maxt);
2117
2118        const Vector3 origin = entp;
2119
2120        VspNode *node = mRoot;
2121        VspNode *farChild = NULL;
2122
2123        float position;
2124        int axis;
2125
2126        while (1)
2127        {
2128                if (!node->IsLeaf())
2129                {
2130                        VspInterior *in = dynamic_cast<VspInterior *>(node);
2131                        position = in->GetPosition();
2132                        axis = in->GetAxis();
2133
2134                        if (entp[axis] <= position)
2135                        {
2136                                if (extp[axis] <= position)
2137                                {
2138                                        node = in->GetBack();
2139                                        // cases N1,N2,N3,P5,Z2,Z3
2140                                        continue;
2141                                }
2142                                else
2143                                {
2144                                        // case N4
2145                                        node = in->GetBack();
2146                                        farChild = in->GetFront();
2147                                }
2148                        }
2149                        else
2150                        {
2151                                if (position <= extp[axis])
2152                                {
2153                                        node = in->GetFront();
2154                                        // cases P1,P2,P3,N5,Z1
2155                                        continue;
2156                                }
2157                                else
2158                                {
2159                                        node = in->GetFront();
2160                                        farChild = in->GetBack();
2161                                        // case P4
2162                                }
2163                        }
2164
2165                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2166                        // case N4 or P4
2167                        const float tdist = (position - origin[axis]) / dir[axis];
2168                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2169                        extp = origin + dir * tdist;
2170                        maxt = tdist;
2171                }
2172                else
2173                {
2174                        // compute intersection with all objects in this leaf
2175                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2176                        ViewCell *vc = leaf->GetViewCell();
2177
2178                        if (!vc->Mailed())
2179                        {
2180                                vc->Mail();
2181                                // todo: add view cells to ray
2182                                ++ hits;
2183                        }
2184#if 0
2185                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2186#endif
2187                        // get the next node from the stack
2188                        if (tStack.empty())
2189                                break;
2190
2191                        entp = extp;
2192                        mint = maxt;
2193                       
2194                        LineTraversalData &s  = tStack.top();
2195                        node = s.mNode;
2196                        extp = s.mExitPoint;
2197                        maxt = s.mMaxT;
2198                        tStack.pop();
2199                }
2200        }
2201
2202        return hits;
2203}
2204
2205
2206ViewCell *VspTree::GetViewCell(const Vector3 &point, const bool active)
2207{
2208        if (mRoot == NULL)
2209                return NULL;
2210
2211        stack<VspNode *> nodeStack;
2212        nodeStack.push(mRoot);
2213 
2214        ViewCellLeaf *viewcell = NULL;
2215 
2216        while (!nodeStack.empty()) 
2217        {
2218                VspNode *node = nodeStack.top();
2219                nodeStack.pop();
2220       
2221                if (node->IsLeaf())
2222                {
2223                        viewcell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
2224                        break;
2225                }
2226                else   
2227                {       
2228                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2229     
2230                        // random decision
2231                        if (interior->GetPosition() - point[interior->GetAxis()] < 0)
2232                                nodeStack.push(interior->GetBack());
2233                        else
2234                                nodeStack.push(interior->GetFront());
2235                }
2236        }
2237 
2238        if (active)
2239                return mViewCellsTree->GetActiveViewCell(viewcell);
2240        else
2241                return viewcell;
2242}
2243
2244
2245bool VspTree::ViewPointValid(const Vector3 &viewPoint) const
2246{
2247        VspNode *node = mRoot;
2248
2249        while (1)
2250        {
2251                // early exit
2252                if (node->TreeValid())
2253                        return true;
2254
2255                if (node->IsLeaf())
2256                        return false;
2257                       
2258                VspInterior *in = dynamic_cast<VspInterior *>(node);
2259                                       
2260                if (in->GetPosition() - viewPoint[in->GetAxis()] <= 0)
2261                {
2262                        node = in->GetBack();
2263                }
2264                else
2265                {
2266                        node = in->GetFront();
2267                }
2268        }
2269
2270        // should never come here
2271        return false;
2272}
2273
2274
2275void VspTree::PropagateUpValidity(VspNode *node)
2276{
2277        const bool isValid = node->TreeValid();
2278
2279        // propagative up invalid flag until only invalid nodes exist over this node
2280        if (!isValid)
2281        {
2282                while (!node->IsRoot() && node->GetParent()->TreeValid())
2283                {
2284                        node = node->GetParent();
2285                        node->SetTreeValid(false);
2286                }
2287        }
2288        else
2289        {
2290                // propagative up valid flag until one of the subtrees is invalid
2291                while (!node->IsRoot() && !node->TreeValid())
2292                {
2293            node = node->GetParent();
2294                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2295                       
2296                        // the parent is valid iff both leaves are valid
2297                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
2298                                                           interior->GetFront()->TreeValid());
2299                }
2300        }
2301}
2302
2303#if ZIPPED_VIEWCELLS
2304bool VspTree::Export(ogzstream &stream)
2305#else
2306bool VspTree::Export(ofstream &stream)
2307#endif
2308{
2309        ExportNode(mRoot, stream);
2310
2311        return true;
2312}
2313
2314
2315#if ZIPPED_VIEWCELLS
2316void VspTree::ExportNode(VspNode *node, ogzstream &stream)
2317#else
2318void VspTree::ExportNode(VspNode *node, ofstream &stream)
2319#endif
2320{
2321        if (node->IsLeaf())
2322        {
2323                VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2324                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
2325
2326                int id = -1;
2327                if (viewCell != mOutOfBoundsCell)
2328                        id = viewCell->GetId();
2329
2330                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
2331        }
2332        else
2333        {       
2334                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2335       
2336                AxisAlignedPlane plane = interior->GetPlane();
2337                stream << "<Interior plane=\"" << plane.mPosition << " "
2338                           << plane.mAxis << "\">" << endl;
2339
2340                ExportNode(interior->GetBack(), stream);
2341                ExportNode(interior->GetFront(), stream);
2342
2343                stream << "</Interior>" << endl;
2344        }
2345}
2346
2347
2348int VspTree::SplitRays(const AxisAlignedPlane &plane,
2349                                           RayInfoContainer &rays,
2350                                           RayInfoContainer &frontRays,
2351                                           RayInfoContainer &backRays) const
2352{
2353        int splits = 0;
2354
2355        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2356
2357        for (rit = rays.begin(); rit != rit_end; ++ rit)
2358        {
2359                RayInfo bRay = *rit;
2360               
2361                VssRay *ray = bRay.mRay;
2362                float t;
2363
2364                // get classification and receive new t
2365                //-- test if start point behind or in front of plane
2366                const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
2367                       
2368
2369#if 1
2370                if (side == 0)
2371                {
2372                        ++ splits;
2373
2374                        if (ray->HasPosDir(plane.mAxis))
2375                        {
2376                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2377                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2378                        }
2379                        else
2380                        {
2381                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2382                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2383                        }
2384                }
2385                else if (side == 1)
2386                {
2387                        frontRays.push_back(bRay);
2388                }
2389                else
2390                {
2391                        backRays.push_back(bRay);
2392                }
2393#else
2394                if (side == 0)
2395                {
2396                        ++ splits;
2397
2398                        if (ray->HasPosDir(plane.mAxis))
2399                        {
2400                                backRays.push_back(RayInfo(ray, bRay.GetMaxT(), t));
2401                                frontRays.push_back(RayInfo(ray, t, bRay.GetMinT()));
2402                        }
2403                        else
2404                        {
2405                                frontRays.push_back(RayInfo(ray, bRay.GetMaxT(), t));
2406                                backRays.push_back(RayInfo(ray, t, bRay.GetMinT()));
2407                        }
2408                }
2409                else if (side == 1)
2410                {
2411                        backRays.push_back(bRay);
2412                }
2413                else
2414                {
2415                        frontRays.push_back(bRay);
2416                       
2417                }
2418#endif
2419        }
2420
2421        return splits;
2422}
2423
2424
2425AxisAlignedBox3 VspTree::GetBBox(VspNode *node) const
2426{
2427        if (!node->GetParent())
2428                return mBoundingBox;
2429
2430        if (!node->IsLeaf())
2431        {
2432                return (dynamic_cast<VspInterior *>(node))->GetBoundingBox();           
2433        }
2434
2435        VspInterior *parent = dynamic_cast<VspInterior *>(node->GetParent());
2436
2437        AxisAlignedBox3 box(parent->GetBoundingBox());
2438
2439        if (parent->GetFront() == node)
2440      box.SetMin(parent->GetAxis(), parent->GetPosition());
2441    else
2442      box.SetMax(parent->GetAxis(), parent->GetPosition());
2443
2444        return box;
2445}
2446
2447
2448
2449
2450int VspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box,
2451                                                                         ViewCellContainer &viewCells) const
2452{
2453        stack<VspNode *> nodeStack;
2454 
2455        ViewCell::NewMail();
2456
2457        while (!nodeStack.empty())
2458        {
2459                VspNode *node = nodeStack.top();
2460                nodeStack.pop();
2461
2462                const AxisAlignedBox3 bbox = GetBBox(node);
2463
2464                if (bbox.Includes(box))
2465                {
2466                        // node geometry is contained in box
2467                        CollectViewCells(node, true, viewCells, true);
2468                }
2469                else if (Overlap(bbox, box))
2470                {
2471                        if (node->IsLeaf())
2472                        {
2473                                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2474                       
2475                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
2476                                {
2477                                        leaf->GetViewCell()->Mail();
2478                                        viewCells.push_back(leaf->GetViewCell());
2479                                }
2480                        }
2481                        else
2482                        {
2483                                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2484                       
2485                                VspNode *first = interior->GetFront();
2486                                VspNode *second = interior->GetBack();
2487           
2488                                nodeStack.push(first);
2489                                nodeStack.push(second);
2490                        }
2491                }       
2492                // default: cull
2493        }
2494
2495        return (int)viewCells.size();
2496}
2497
2498
2499
2500/*****************************************************************/
2501/*                class OspTree implementation                   */
2502/*****************************************************************/
2503
2504
2505OspTree::OspTree():
2506mRoot(NULL)
2507#if TODO
2508mOutOfBoundsCell(NULL),
2509mStoreRays(false),
2510mUseRandomAxis(false),
2511mTimeStamp(1)
2512#endif
2513{
2514#if TODO
2515        bool randomize = false;
2516        Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
2517        if (randomize)
2518                Randomize(); // initialise random generator for heuristics
2519
2520        //-- termination criteria for autopartition
2521        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxDepth", mTermMaxDepth);
2522        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minPvs", mTermMinPvs);
2523        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minRays", mTermMinRays);
2524        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minProbability", mTermMinProbability);
2525        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxRayContribution", mTermMaxRayContribution);
2526       
2527        Environment::GetSingleton()->GetIntValue("VspTree.Termination.missTolerance", mTermMissTolerance);
2528        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxViewCells", mMaxViewCells);
2529
2530        //-- max cost ratio for early tree termination
2531        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxCostRatio", mTermMaxCostRatio);
2532
2533        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
2534        Environment::GetSingleton()->GetIntValue("VspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
2535
2536        // HACK//mTermMinPolygons = 25;
2537
2538        //-- factors for bsp tree split plane heuristics
2539        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.ct_div_ci", mCtDivCi);
2540
2541        //-- partition criteria
2542        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.epsilon", mEpsilon);
2543
2544        // if only the driving axis is used for axis aligned split
2545        Environment::GetSingleton()->GetBoolValue("VspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
2546       
2547        //Environment::GetSingleton()->GetFloatValue("VspTree.maxTotalMemory", mMaxTotalMemory);
2548        Environment::GetSingleton()->GetFloatValue("VspTree.maxStaticMemory", mMaxMemory);
2549
2550        Environment::GetSingleton()->GetBoolValue("VspTree.useCostHeuristics", mUseCostHeuristics);
2551        Environment::GetSingleton()->GetBoolValue("VspTree.simulateOctree", mCirculatingAxis);
2552
2553
2554        char subdivisionStatsLog[100];
2555        Environment::GetSingleton()->GetStringValue("VspTree.subdivisionStats", subdivisionStatsLog);
2556        mSubdivisionStats.open(subdivisionStatsLog);
2557
2558        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.minBand", mMinBand);
2559        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.maxBand", mMaxBand);
2560       
2561
2562        //-- debug output
2563
2564        Debug << "******* OSP options ******** " << endl;
2565    Debug << "max depth: " << mTermMaxDepth << endl;
2566        Debug << "min PVS: " << mTermMinPvs << endl;
2567        Debug << "min probabiliy: " << mTermMinProbability << endl;
2568        Debug << "min rays: " << mTermMinRays << endl;
2569        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
2570        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
2571        Debug << "miss tolerance: " << mTermMissTolerance << endl;
2572        Debug << "max view cells: " << mMaxViewCells << endl;
2573        Debug << "randomize: " << randomize << endl;
2574
2575        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
2576        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
2577        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
2578        Debug << "max memory: " << mMaxMemory << endl;
2579        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
2580        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
2581       
2582        Debug << "circulating axis: " << mCirculatingAxis << endl;
2583        Debug << "minband: " << mMinBand << endl;
2584        Debug << "maxband: " << mMaxBand << endl;
2585       
2586
2587        mSplitCandidates = new vector<SortableEntry>;
2588
2589        Debug << endl;
2590#endif
2591}
2592
2593
2594
2595void OspTree::SplitObjects(const AxisAlignedPlane & splitPlane,
2596                                                   const ObjectContainer &objects,
2597                                                   ObjectContainer &back,
2598                                                   ObjectContainer &front)
2599{
2600        ObjectContainer::const_iterator oit, oit_end = objects.end();
2601
2602    for (oit = objects.begin(); oit != oit_end; ++ oit)
2603        {
2604                // determine the side of this ray with respect to the plane
2605                const AxisAlignedBox3 box = (*oit)->GetBox();
2606
2607                if (box.Max(splitPlane.mAxis) >= splitPlane.mPosition)
2608            front.push_back(*oit);
2609   
2610                if (box.Min(splitPlane.mAxis) < splitPlane.mPosition)
2611                        back.push_back(*oit);
2612#if TODO   
2613                mStat.objectRefs -= (int)objects.size();
2614                mStat.objectRefs += objectsBack + objectsFront;
2615#endif
2616        }
2617}
2618
2619
2620KdInterior *OspTree::SubdivideNode(KdLeaf *leaf,
2621                                                                   const AxisAlignedPlane &splitPlane,
2622                                                                   const AxisAlignedBox3 &box,
2623                                                                   AxisAlignedBox3 &backBBox,
2624                                                                   AxisAlignedBox3 &frontBBox)
2625{
2626#if TODO
2627    mSpatialStat.nodes += 2;
2628        mSpatialStat.splits[axis];
2629#endif
2630
2631        // add the new nodes to the tree
2632        KdInterior *node = new KdInterior(leaf->mParent);
2633
2634        const int axis = splitPlane.mAxis;
2635        const float position = splitPlane.mPosition;
2636
2637        node->mAxis = axis;
2638        node->mPosition = position;
2639        node->mBox = box;
2640
2641    backBBox = box;
2642        frontBBox = box;
2643 
2644        // first count ray sides
2645        int objectsBack = 0;
2646        int objectsFront = 0;
2647 
2648        backBBox.SetMax(axis, position);
2649        frontBBox.SetMin(axis, position);
2650
2651        ObjectContainer::const_iterator mi, mi_end = leaf->mObjects.end();
2652
2653        for ( mi = leaf->mObjects.begin(); mi != mi_end; ++ mi)
2654        {
2655                // determine the side of this ray with respect to the plane
2656                const AxisAlignedBox3 box = (*mi)->GetBox();
2657               
2658                if (box.Max(axis) > position)
2659                        ++ objectsFront;
2660   
2661                if (box.Min(axis) < position)
2662                        ++ objectsBack;
2663        }
2664
2665        KdLeaf *back = new KdLeaf(node, objectsBack);
2666        KdLeaf *front = new KdLeaf(node, objectsFront);
2667
2668        // replace a link from node's parent
2669        if (leaf->mParent)
2670                leaf->mParent->ReplaceChildLink(leaf, node);
2671
2672        // and setup child links
2673        node->SetupChildLinks(back, front);
2674
2675        SplitObjects(splitPlane, leaf->mObjects, back->mObjects, front->mObjects);
2676 
2677        //delete leaf;
2678        return node;
2679}
2680
2681
2682KdNode *OspTree::Subdivide(SplitQueue &tQueue,
2683                                                   OspSplitCandidate &splitCandidate,
2684                                                   const bool globalCriteriaMet)
2685{
2686        OspTraversalData &tData = splitCandidate.mParentData;
2687
2688        KdNode *newNode = tData.mNode;
2689
2690        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
2691        {       
2692                OspTraversalData tFrontData;
2693                OspTraversalData tBackData;
2694
2695                //-- continue subdivision
2696               
2697                // create new interior node and two leaf node
2698                const AxisAlignedPlane splitPlane = splitCandidate.mSplitPlane;
2699               
2700                newNode = SubdivideNode(dynamic_cast<KdLeaf *>(newNode),
2701                                                                splitPlane,
2702                                                                tData.mBoundingBox,
2703                                                                tFrontData.mBoundingBox,
2704                                                                tBackData.mBoundingBox);
2705       
2706                const int maxCostMisses = splitCandidate.mMaxCostMisses;
2707
2708                // how often was max cost ratio missed in this branch?
2709                tFrontData.mMaxCostMisses = maxCostMisses;
2710                tBackData.mMaxCostMisses = maxCostMisses;
2711                       
2712                //-- push the new split candidates on the queue
2713                OspSplitCandidate *frontCandidate = new OspSplitCandidate();
2714                OspSplitCandidate *backCandidate = new OspSplitCandidate();
2715
2716                EvalSplitCandidate(tFrontData, *frontCandidate);
2717                EvalSplitCandidate(tBackData, *backCandidate);
2718
2719                tQueue.push(frontCandidate);
2720                tQueue.push(backCandidate);
2721
2722                // delete old leaf node
2723                DEL_PTR(tData.mNode);
2724        }
2725
2726
2727        //-- terminate traversal
2728        if (newNode->IsLeaf())
2729        {
2730                //KdLeaf *leaf = dynamic_cast<KdLeaf *>(newNode);
2731                EvaluateLeafStats(tData);               
2732        }
2733
2734        //-- cleanup
2735        tData.Clear();
2736       
2737        return newNode;
2738}
2739
2740
2741void OspTree::EvalSplitCandidate(OspTraversalData &tData,
2742                                                                 OspSplitCandidate &splitData)
2743{
2744        float frontProb;
2745        float backtProb;
2746       
2747        KdLeaf *leaf = dynamic_cast<KdLeaf *>(tData.mNode);
2748
2749        // compute locally best split plane
2750        const bool success = false;
2751#if TODO
2752        SelectPlane(tData, splitData.mSplitPlane,
2753                                frontData.mProbability, backData.mProbability);
2754
2755        //TODO
2756        // compute global decrease in render cost
2757        splitData.mPriority = EvalRenderCostDecrease(splitData.mSplitPlane, tData);
2758        splitData.mParentData = tData;
2759        splitData.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
2760#endif
2761}
2762
2763
2764bool OspTree::LocalTerminationCriteriaMet(const OspTraversalData &data) const
2765{
2766        // matt: TODO
2767        return true;
2768    /*          (((int)data.mRays->size() <= mTermMinRays) ||
2769                 (data.mPvs <= mTermMinPvs)   ||
2770                 (data.mProbability <= mTermMinProbability) ||
2771                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
2772                 (data.mDepth >= mTermMaxDepth));*/
2773}
2774
2775
2776bool OspTree::GlobalTerminationCriteriaMet(const OspTraversalData &data) const
2777{
2778        // matt: TODO
2779        return true;
2780                /*(mOutOfMemory ||
2781                (mVspStats.Leaves() >= mMaxViewCells) ||
2782                (mGlobalCostMisses >= mTermGlobalCostMissTolerance));*/
2783}
2784
2785
2786void OspTree::EvaluateLeafStats(const OspTraversalData &data)
2787{
2788#if TODO
2789        // the node became a leaf -> evaluate stats for leafs
2790        VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);
2791
2792        if (data.mPvs > mVspStats.maxPvs)
2793        {
2794                mVspStats.maxPvs = data.mPvs;
2795        }
2796
2797        mVspStats.pvs += data.mPvs;
2798
2799        if (data.mDepth < mVspStats.minDepth)
2800        {
2801                mVspStats.minDepth = data.mDepth;
2802        }
2803       
2804        if (data.mDepth >= mTermMaxDepth)
2805        {
2806        ++ mVspStats.maxDepthNodes;
2807                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
2808        }
2809
2810        // accumulate rays to compute rays /  leaf
2811        mVspStats.accumRays += (int)data.mRays->size();
2812
2813        if (data.mPvs < mTermMinPvs)
2814                ++ mVspStats.minPvsNodes;
2815
2816        if ((int)data.mRays->size() < mTermMinRays)
2817                ++ mVspStats.minRaysNodes;
2818
2819        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
2820                ++ mVspStats.maxRayContribNodes;
2821
2822        if (data.mProbability <= mTermMinProbability)
2823                ++ mVspStats.minProbabilityNodes;
2824       
2825        // accumulate depth to compute average depth
2826        mVspStats.accumDepth += data.mDepth;
2827
2828        ++ mCreatedViewCells;
2829
2830#ifdef _DEBUG
2831        Debug << "BSP stats: "
2832                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
2833                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
2834                  << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
2835                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
2836                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "), "
2837                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
2838#endif
2839
2840#endif
2841}
2842
2843
2844/*
2845float OspTree::EvalLocalCostHeuristics(const RayInfoContainer &rays,
2846                                                                           const AxisAlignedBox3 &box,
2847                                                                           int pvsSize,
2848                                                                           const int axis,
2849                                                                           float &position)
2850{
2851        const float minBox = box.Min(axis);
2852        const float maxBox = box.Max(axis);
2853
2854        const float sizeBox = maxBox - minBox;
2855
2856        const float minBand = minBox + mMinBand * sizeBox;
2857        const float maxBand = minBox + mMaxBand * sizeBox;
2858
2859        SortSplitCandidates(rays, axis, minBand, maxBand);
2860
2861        // prepare the sweep
2862        // (note: returns pvs size, so there would be no need
2863        // to give pvs size as argument)
2864        pvsSize = PrepareHeuristics(rays);
2865
2866        Debug << "here45 pvs: " << pvsSize << endl;
2867
2868        // go through the lists, count the number of objects left and right
2869        // and evaluate the following cost funcion:
2870        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
2871
2872        int pvsl = 0;
2873        int pvsr = pvsSize;
2874
2875        int pvsBack = pvsl;
2876        int pvsFront = pvsr;
2877
2878        float sum = (float)pvsSize * sizeBox;
2879        float minSum = 1e20f;
2880
2881       
2882        // if no good split can be found, take mid split
2883        position = minBox + 0.5f * sizeBox;
2884       
2885        // the relative cost ratio
2886        float ratio = 99999999.0f;
2887        bool splitPlaneFound = false;
2888
2889        Intersectable::NewMail();
2890        KdLeaf::NewMail();
2891
2892
2893        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
2894
2895        //-- traverse through visibility events
2896        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
2897        {
2898                EvalPvsIncr(*ci, pvsl, pvsr);
2899
2900                // Note: sufficient to compare size of bounding boxes of front and back side?
2901                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
2902                {
2903                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
2904
2905                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << "\t cost= " << sum << endl;
2906
2907                        if (sum < minSum)
2908                        {
2909                                splitPlaneFound = true;
2910
2911                                minSum = sum;
2912                                position = (*ci).value;
2913                               
2914                                pvsBack = pvsl;
2915                                pvsFront = pvsr;
2916                        }
2917                }
2918        }
2919       
2920       
2921        // -- compute cost
2922
2923        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
2924        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
2925
2926        const float pOverall = sizeBox;
2927        const float pBack = position - minBox;
2928        const float pFront = maxBox - position;
2929
2930        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
2931    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
2932        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
2933       
2934        const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
2935        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
2936
2937        if (splitPlaneFound)
2938        {
2939                ratio = newRenderCost / oldRenderCost;
2940        }
2941       
2942        //if (axis != 1)
2943        Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
2944              <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
2945
2946        return ratio;
2947}
2948*/
2949
2950float OspTree::EvalLocalCostHeuristics(BspLeaf *node,
2951                                                                           const AxisAlignedBox3 &box,
2952                                                                           const int axis,
2953                                                                           float &position,
2954                                                                           int &objectsBack,
2955                                                                           int &objectsFront)
2956{
2957       
2958        /*SortSplitCandidates(node, axis);
2959 
2960        // go through the lists, count the number of objects left and right
2961        // and evaluate the following cost funcion:
2962        // C = ct_div_ci  + (ol + or)/queries
2963
2964        float totalIntersections = 0.0f;
2965       
2966        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
2967
2968        for(ci = splitCandidates->begin(); ci != ci_end; ++ ci)
2969        {
2970                 if ((*ci).type == SortableEntry::BOX_MIN)
2971                 {
2972                         totalIntersections += (*ci).intersectable->IntersectionComplexity();
2973                 }
2974        }
2975       
2976        float intersectionsLeft = 0;
2977        float intersectionsRight = totalIntersections;
2978
2979       
2980        int objectsLeft = 0, objectsRight = (int)node->mObjects.size();
2981 
2982        float minBox = box.Min(axis);
2983        float maxBox = box.Max(axis);
2984       
2985        float boxArea = box.SurfaceArea();
2986 
2987        float minBand = minBox + mSplitBorder*(maxBox - minBox);
2988        float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
2989
2990    float minSum = 1e20f;
2991 
2992        for(ci = splitCandidates->begin(); ci != ci_end; ++ ci)
2993        {
2994                EvalPvsIncr(ci, pvsl, pvsr);
2995
2996                if ((*ci).value > minBand && (*ci).value < maxBand)
2997                {
2998                        AxisAlignedBox3 lbox = box;
2999                        AxisAlignedBox3 rbox = box;
3000                       
3001                        lbox.SetMax(axis, (*ci).value);
3002                        rbox.SetMin(axis, (*ci).value);
3003
3004                        float sum;
3005                       
3006                        if (mSahUseFaces)
3007                                sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
3008                        else
3009                                sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
3010     
3011                        //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
3012                        //      cout<<"cost= "<<sum<<endl;
3013     
3014                        if (sum < minSum)
3015                        {
3016                                minSum = sum;
3017                                position = (*ci).value;
3018
3019                                objectsBack = objectsLeft;
3020                                objectsFront = objectsRight;
3021                        }
3022                }
3023        }
3024 
3025        float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
3026        float newCost = mCt_div_ci + minSum/boxArea;
3027        float ratio = newCost/oldCost;
3028 
3029#if 0
3030        cout<<"===================="<<endl;
3031        cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
3032                <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
3033#endif
3034       
3035        return ratio;
3036        */
3037        return 0;
3038}
3039
3040
3041void OspTree::SortSplitCandidates(const RayInfoContainer &rays,
3042                                                                  const int axis,
3043                                                                  float minBand,
3044                                                                  float maxBand)
3045{
3046        // TODO
3047}
3048
3049
3050
3051int OspTree::PrepareHeuristics(Intersectable *object)
3052{
3053        ViewCellPvsMap::const_iterator vit, vit_end = object->mViewCellPvs.mEntries.end();
3054       
3055        int pvsSize = 0;
3056
3057        for (vit = object->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3058        {
3059        ViewCell *vc = (*vit).first;
3060
3061                if (!vc->Mailed())
3062                {
3063                        vc->Mail();
3064                        vc->mCounter = 1;
3065                        ++ pvsSize;
3066                }
3067                else
3068                {
3069                        ++ vc->mCounter;
3070                }
3071        }
3072
3073        return pvsSize;
3074}
3075
3076
3077int OspTree::PrepareHeuristics(const RayInfoContainer &rays)
3078{       
3079        Intersectable::NewMail();
3080        KdNode::NewMail();
3081
3082        int pvsSize = 0;
3083
3084        RayInfoContainer::const_iterator ri, ri_end = rays.end();
3085
3086        //-- set all kd nodes as belonging to the front pvs
3087
3088        for (ri = rays.begin(); ri != ri_end; ++ ri)
3089        {
3090                Intersectable *oObject = (*ri).mRay->mOriginObject;
3091
3092                if (oObject)
3093                {
3094                        pvsSize += PrepareHeuristics(oObject);         
3095                }
3096
3097                Intersectable *tObject = (*ri).mRay->mTerminationObject;
3098
3099                if (tObject)
3100                {
3101                        pvsSize += PrepareHeuristics(tObject); 
3102                }
3103        }
3104
3105        return pvsSize;
3106}
3107
3108
3109void OspTree::EvalPvsIncr(const SortableEntry &ci,
3110                                                  int &pvsLeft,
3111                                                  int &pvsRight) const
3112{
3113        Intersectable *obj = ci.mObject;
3114
3115        switch (ci.type)
3116        {
3117                case SortableEntry::BOX_MIN:
3118                        AddContriToPvs(obj, pvsLeft);
3119                        break;
3120                       
3121                case SortableEntry::BOX_MAX:
3122                        RemoveContriFromPvs(obj, pvsRight);
3123                        break;
3124        }
3125}
3126
3127
3128void OspTree::RemoveContriFromPvs(Intersectable *object, int &pvs) const
3129{
3130        ViewCellPvsMap::const_iterator vit, vit_end = object->mViewCellPvs.mEntries.end();
3131       
3132        for (vit = object->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3133        {
3134        ViewCell *vc = (*vit).first;
3135
3136                if (-- vc->mCounter == 0)
3137                {
3138                        -- pvs;
3139                }
3140        }
3141}
3142
3143
3144void OspTree::AddContriToPvs(Intersectable *object, int &pvs) const
3145{
3146        ViewCellPvsMap::const_iterator vit, vit_end = object->mViewCellPvs.mEntries.end();
3147
3148        for (vit = object->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3149        {
3150        ViewCell *vc = (*vit).first;
3151
3152                if (!vc->Mailed())
3153                {
3154                        vc->Mail();             
3155                        ++ pvs;
3156                }
3157        }
3158}
3159
3160
3161/*********************************************************************/
3162/*                class HierarchyManager implementation              */
3163/*********************************************************************/
3164
3165
3166
3167HierarchyManager::HierarchyManager(VspTree &vspTree, OspTree &ospTree):
3168mVspTree(vspTree), mOspTree(ospTree)
3169{
3170}
3171
3172
3173SplitCandidate *HierarchyManager::NextSplitCandidate()
3174{
3175        SplitCandidate *splitCandidate = mTQueue.top();
3176        //Debug << "priority: " << splitCandidate->GetPriority() << endl;
3177        mTQueue.pop();
3178
3179        return splitCandidate;
3180}
3181
3182
3183void HierarchyManager::PrepareConstruction(const VssRayContainer &sampleRays,
3184                                                                                   const ObjectContainer &objects,
3185                                                                                   AxisAlignedBox3 *forcedViewSpace,
3186                                                                                   RayInfoContainer &rays)
3187{
3188        mVspTree.PrepareConstruction(sampleRays, forcedViewSpace);
3189
3190        long startTime = GetTime();
3191
3192        cout << "storing rays ... ";
3193
3194        Intersectable::NewMail();
3195
3196        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
3197
3198        //-- store rays
3199        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
3200        {
3201                VssRay *ray = *rit;
3202
3203                float minT, maxT;
3204
3205                static Ray hray;
3206                hray.Init(*ray);
3207               
3208                // TODO: not very efficient to implictly cast between rays types
3209                if (mVspTree.GetBoundingBox().GetRaySegment(hray, minT, maxT))
3210                {
3211                        float len = ray->Length();
3212
3213                        if (!len)
3214                                len = Limits::Small;
3215
3216                        rays.push_back(RayInfo(ray, minT / len, maxT / len));
3217                }
3218        }
3219
3220        cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
3221
3222        /// add first candidate for view space partition
3223        mVspTree.mRoot = new VspLeaf();
3224        const float prop = mVspTree.mBoundingBox.GetVolume();
3225
3226        /// first traversal data
3227        VspTree::VspTraversalData tData(mVspTree.mRoot,
3228                                                                        0,
3229                                                                        &rays,
3230                                                                        //(int)objects.size(),
3231                                                                        mVspTree.ComputePvsSize(rays),
3232                                                                        prop,
3233                                                                        mVspTree.mBoundingBox);
3234
3235        // compute first split candidate
3236        VspTree::VspSplitCandidate *splitCandidate = new VspTree::VspSplitCandidate();
3237    mVspTree.EvalSplitCandidate(tData, *splitCandidate);
3238
3239        mTQueue.push(splitCandidate);
3240
3241        //mOspTree->PrepareConstruction(sampleRays, forcedViewSpace, rays);
3242}
3243
3244
3245bool HierarchyManager::GlobalTerminationCriteriaMet(SplitCandidate *candidate) const
3246{
3247        if (candidate->Type() == SplitCandidate::VIEW_SPACE)
3248        {
3249                VspTree::VspSplitCandidate *sc =
3250                        dynamic_cast<VspTree::VspSplitCandidate *>(candidate);
3251               
3252                return mVspTree.GlobalTerminationCriteriaMet(sc->mParentData);
3253        }
3254
3255        return true;
3256}
3257
3258
3259void HierarchyManager::Construct(const VssRayContainer &sampleRays,
3260                                                                 const ObjectContainer &objects,
3261                                                                 AxisAlignedBox3 *forcedViewSpace)
3262{
3263        RayInfoContainer *rays = new RayInfoContainer();
3264       
3265        // prepare vsp and osp trees for traversal
3266        PrepareConstruction(sampleRays, objects, forcedViewSpace, *rays);
3267
3268        mVspTree.mVspStats.Reset();
3269        mVspTree.mVspStats.Start();
3270
3271        cout << "Constructing view space / object space tree ... \n";
3272        const long startTime = GetTime();       
3273       
3274        int i = 0;
3275        while (!FinishedConstruction())
3276        {
3277                SplitCandidate *splitCandidate = NextSplitCandidate();
3278           
3279                const bool globalTerminationCriteriaMet =
3280                        GlobalTerminationCriteriaMet(splitCandidate);
3281
3282                cout << "view cells: " << i ++ << endl;
3283
3284                // cost ratio of cost decrease / totalCost
3285                const float costRatio = splitCandidate->GetPriority() / mTotalCost;
3286                //Debug << "cost ratio: " << costRatio << endl;
3287
3288                if (costRatio < mTermMinGlobalCostRatio)
3289                        ++ mGlobalCostMisses;
3290
3291                //-- subdivide leaf node
3292                //-- either a object space or view space split
3293                if (splitCandidate->Type() == SplitCandidate::VIEW_SPACE)
3294                {
3295                        VspTree::VspSplitCandidate *sc =
3296                                dynamic_cast<VspTree::VspSplitCandidate *>(splitCandidate);
3297
3298                        VspNode *r = mVspTree.Subdivide(mTQueue, *sc, globalTerminationCriteriaMet);
3299                }
3300                else // object space split
3301                {
3302#if TODO
3303                        KdNode *r = mKdtree.Subdivide(tOspQueue, dynamic_cast<OspSplitCandidate<(splitCandidate));
3304#endif
3305                }
3306
3307                DEL_PTR(splitCandidate);
3308        }
3309
3310        cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl;
3311
3312        mVspTree.mVspStats.Stop();
3313}
3314
3315
3316bool HierarchyManager::FinishedConstruction()
3317{
3318        return mTQueue.empty();
3319}
3320
3321
3322}
Note: See TracBrowser for help on using the repository browser.