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

Revision 1093, 85.2 KB checked in by mattausch, 19 years ago (diff)
Line 
1#include <stack>
2#include <time.h>
3#include <iomanip>
4
5#include "ViewCell.h"
6#include "Plane3.h"
7#include "VspOspTree.h"
8#include "Mesh.h"
9#include "common.h"
10#include "Environment.h"
11#include "Polygon3.h"
12#include "Ray.h"
13#include "AxisAlignedBox3.h"
14#include "Exporter.h"
15#include "Plane3.h"
16#include "ViewCellsManager.h"
17#include "Beam.h"
18#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        interior->mTimeStamp = mTimeStamp ++;
732       
733        return interior;
734}
735
736
737
738void VspTree::ProcessViewCellObjects(ViewCell *parent, ViewCell *front, ViewCell *back) const
739{
740        if (parent)
741        {
742                // remove the parents from the object pvss
743                ObjectPvsMap::const_iterator oit, oit_end = parent->GetPvs().mEntries.end();
744
745                for (oit = parent->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
746                {
747                        Intersectable *object = (*oit).first;
748                        // HACK: make sure that the view cell is removed from the pvs
749                        const float high_contri = 99999999999;
750                        object->mViewCellPvs.RemoveSample(parent, 999999);
751                }
752        }
753
754        if (front)
755        {
756                // Add front view cell to the object pvsss
757                ObjectPvsMap::const_iterator oit, oit_end = front->GetPvs().mEntries.end();
758
759                for (oit = front->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
760                {
761                        Intersectable *object = (*oit).first;
762                        object->mViewCellPvs.AddSample(front, 1);
763                }
764        }
765
766        if (back)
767        {
768                // Add back view cell to the object pvsss
769                ObjectPvsMap::const_iterator oit, oit_end = back->GetPvs().mEntries.end();
770
771                for (oit = back->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
772                {
773                        Intersectable *object = (*oit).first;
774                        object->mViewCellPvs.AddSample(back, 1);
775                }
776        }
777}
778
779
780void VspTree::AddToPvs(VspLeaf *leaf,
781                                           const RayInfoContainer &rays,
782                                           float &sampleContributions,
783                                           int &contributingSamples)
784{
785        sampleContributions = 0;
786        contributingSamples = 0;
787 
788        RayInfoContainer::const_iterator it, it_end = rays.end();
789 
790        ViewCellLeaf *vc = leaf->GetViewCell();
791 
792        // add contributions from samples to the PVS
793        for (it = rays.begin(); it != it_end; ++ it)
794        {
795                float sc = 0.0f;
796                VssRay *ray = (*it).mRay;
797
798                bool madeContrib = false;
799                float contribution;
800
801                if (ray->mTerminationObject)
802                {
803                        if (vc->AddPvsSample(ray->mTerminationObject, ray->mPdf, contribution))
804                        {
805                                madeContrib = true;
806                        }
807
808                        sc += contribution;
809                }
810         
811                if (ray->mOriginObject)
812                {
813                        if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
814                        {
815                                madeContrib = true;
816                        }
817
818                        sc += contribution;
819                }
820         
821                sampleContributions += sc;
822
823                if (madeContrib)
824                        ++ contributingSamples;
825               
826                // store rays for visualization
827                if (0) leaf->mVssRays.push_back(new VssRay(*ray));
828        }
829}
830
831
832void VspTree::SortSplitCandidates(const RayInfoContainer &rays,
833                                                                  const int axis,
834                                                                  float minBand,
835                                                                  float maxBand)
836{
837        mSplitCandidates->clear();
838
839        int requestedSize = 2 * (int)(rays.size());
840
841        // creates a sorted split candidates array
842        if (mSplitCandidates->capacity() > 500000 &&
843                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
844        {
845        delete mSplitCandidates;
846                mSplitCandidates = new vector<SortableEntry>;
847        }
848
849        mSplitCandidates->reserve(requestedSize);
850
851        float pos;
852
853        //-- insert all queries
854        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
855        {
856                const bool positive = (*ri).mRay->HasPosDir(axis);
857                               
858                pos = (*ri).ExtrapOrigin(axis);
859               
860                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
861                                                                        pos, (*ri).mRay));
862
863                pos = (*ri).ExtrapTermination(axis);
864
865                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
866                                                                        pos, (*ri).mRay));
867        }
868
869        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
870}
871
872
873int VspTree::GetPvsContribution(Intersectable *object) const
874{
875    int pvsContri = 0;
876
877        KdPvsMap::const_iterator kit, kit_end = object->mKdPvs.mEntries.end();
878
879        Intersectable::NewMail();
880
881        // Search kd leaves this object is attached to
882        for (kit = object->mKdPvs.mEntries.begin(); kit != kit_end; ++ kit)
883        {
884                KdNode *l = (*kit).first;
885
886                // new object found during sweep
887                // => increase pvs contribution of this kd node
888                if (!l->Mailed())
889                {
890                        l->Mail();
891                        ++ pvsContri;
892                }
893        }
894
895        return pvsContri;
896}
897
898
899int VspTree::PrepareHeuristics(Intersectable *object)
900{
901        set<KdLeaf *>::const_iterator kit, kit_end = object->mKdLeaves.end();
902       
903        int pvsSize = 0;
904
905        for (kit = object->mKdLeaves.begin(); kit != kit_end; ++ kit)
906        {
907                KdLeaf *node = dynamic_cast<KdLeaf *>(*kit);
908
909                if (!node->Mailed())
910                {
911                        node->Mail();
912                        node->mCounter = 1;
913
914                        //Debug << "here5 "<<(int)node->mObjects.size() <<" "<< node->mMultipleObjects.size()<< " " <<(int)node->mObjects.size() - node->mMultipleObjects.size()<<endl;
915                        // add objects without the objects which are in several kd leaves
916                        pvsSize += (int)(node->mObjects.size() - node->mMultipleObjects.size());
917                }
918                else
919                {
920                        ++ node->mCounter;
921                }
922
923                //-- the objects belonging to several leaves must be handled seperately
924                ObjectContainer::const_iterator oit, oit_end = node->mMultipleObjects.end();
925
926                for (oit = node->mMultipleObjects.begin(); oit != oit_end; ++ oit)
927                {
928                        Intersectable *object = *oit;
929                                               
930                        if (!object->Mailed())
931                        {
932                                //Debug << "here233: " << object->mKdLeaves.size() << endl;
933                                object->Mail();
934                                object->mCounter = 1;
935
936                                ++ pvsSize;
937                        }
938                        else
939                        {
940                                ++ object->mCounter;
941                        }
942                }
943        }
944
945        return pvsSize;
946}
947
948
949int VspTree::PrepareHeuristics(const RayInfoContainer &rays)
950{       
951        Intersectable::NewMail();
952        KdNode::NewMail();
953
954        int pvsSize = 0;
955
956        RayInfoContainer::const_iterator ri, ri_end = rays.end();
957
958        //-- set all kd nodes as belonging to the front pvs
959
960        for (ri = rays.begin(); ri != ri_end; ++ ri)
961        {
962                Intersectable *oObject = (*ri).mRay->mOriginObject;
963
964                if (oObject)
965                {
966                        if (mPvsCountMethod == PER_OBJECT)
967                        {
968                                if (!oObject->Mailed())
969                                {
970                                        oObject->Mail();
971                                        oObject->mCounter = 1;
972                                        ++ pvsSize;
973                                }
974                                else
975                                {
976                                        ++ oObject->mCounter;
977                                }
978                        }
979                        else
980                        {
981                                pvsSize += PrepareHeuristics(oObject); 
982                        }       
983                }
984
985                Intersectable *tObject = (*ri).mRay->mTerminationObject;
986
987                if (tObject)
988                {
989                        if (mPvsCountMethod == PER_OBJECT)
990                        {
991                                if (!tObject->Mailed())
992                                {
993                                        tObject->Mail();
994                                        tObject->mCounter = 1;
995                                        ++ pvsSize;
996                                }
997                                else
998                                {
999                                        ++ tObject->mCounter;
1000                                }
1001                        }
1002                        else
1003                        {
1004                                pvsSize += PrepareHeuristics(tObject); 
1005                        }       
1006                }
1007        }
1008
1009        return pvsSize;
1010}
1011
1012
1013int VspTree::GetPvsIncr(Intersectable *object, const KdPvsMap &activeNodes)
1014{
1015        // TODO:
1016        // use kd info table to apply set theory in order to
1017        // sort out dublicate pvs entries due to objects which
1018        // belong to more than one kd leaves.
1019#if TODO
1020        KdPvsMap::const_iterator kit, kit_end = obj->mKdPvs.mEntries.end();
1021
1022        // Search kd leaves this object is attached to
1023        for (kit = obj->mKdPvs.mEntries.begin(); kit != kit_end; ++ kit)
1024        {
1025                KdNode *l = (*kit).first;
1026
1027                // new object found during sweep
1028                // => increase pvs contribution of this kd node
1029                if (!l->Mailed())
1030                {
1031                        l->Mail();
1032                        ++ pvsContr;
1033                }
1034        }
1035#endif
1036        return 0;
1037}
1038
1039
1040void VspTree::RemoveContriFromPvs(KdLeaf *leaf, int &pvs) const
1041{
1042        // leaf falls out of right pvs
1043        if (-- leaf->mCounter == 0)
1044        {
1045                pvs -= ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
1046        }
1047
1048        //-- handle objects which are in several kd leaves separately
1049        ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1050
1051        for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1052        {
1053                Intersectable *object = *oit;
1054
1055                if (-- object->mCounter == 0)
1056                {
1057                        -- pvs;
1058                }
1059        }
1060}
1061
1062
1063void VspTree::AddContriToPvs(KdLeaf *leaf, int &pvs) const
1064{
1065        if (!leaf->Mailed())
1066        {
1067                leaf->Mail();
1068
1069                // add objects without those which are part of several kd leaves
1070                pvs += ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
1071
1072                //-- handle objects of several kd leaves separately
1073                ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1074
1075                for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1076                {
1077                        Intersectable *object = *oit;
1078
1079                        // object not previously in left pvs
1080                        if (!object->Mailed())
1081                        {
1082                                object->Mail();
1083                                ++ pvs;
1084                        }
1085                }
1086        }                                       
1087}
1088
1089
1090void VspTree::EvalPvsIncr(const SortableEntry &ci,
1091                                                  int &pvsLeft,
1092                                                  int &pvsRight) const
1093{
1094        VssRay *ray = ci.ray;
1095
1096        Intersectable *oObject = ray->mOriginObject;
1097        Intersectable *tObject = ray->mTerminationObject;
1098
1099        if (oObject)
1100        {       
1101                if (mPvsCountMethod == PER_OBJECT)
1102                {
1103                        if (ci.type == SortableEntry::ERayMin)
1104                        {
1105                                if (!oObject->Mailed())
1106                                {
1107                                        oObject->Mail();
1108                                        ++ pvsLeft;
1109                                }
1110                        }
1111                        else if (ci.type == SortableEntry::ERayMax)
1112                        {
1113                                if (-- oObject->mCounter == 0)
1114                                        -- pvsRight;
1115                        }
1116                }
1117                else // per kd node
1118                {
1119                        set<KdLeaf *>::const_iterator kit, kit_end = oObject->mKdLeaves.end();
1120
1121                        for (kit = oObject->mKdLeaves.begin(); kit != kit_end; ++ kit)
1122                        {
1123                                KdLeaf *node = dynamic_cast<KdLeaf *>(*kit);
1124
1125                                // add contributions of the kd nodes
1126                                if (ci.type == SortableEntry::ERayMin)
1127                                {
1128                                        AddContriToPvs(node, pvsLeft);
1129                                }
1130                                else if (ci.type == SortableEntry::ERayMax)
1131                                {
1132                                        RemoveContriFromPvs(node, pvsRight);
1133                                }
1134                        }
1135                }
1136                       
1137
1138        }
1139       
1140        if (tObject)
1141        {       
1142                if (mPvsCountMethod == PER_OBJECT)
1143                {
1144                        if (ci.type == SortableEntry::ERayMin)
1145                        {
1146                                if (!tObject->Mailed())
1147                                {
1148                                        tObject->Mail();
1149                                        ++ pvsLeft;
1150                                }
1151                        }
1152                        else if (ci.type == SortableEntry::ERayMax)
1153                        {
1154                                if (-- tObject->mCounter == 0)
1155                                        -- pvsRight;
1156                        }
1157                }
1158                else // per kd node
1159                {
1160                        set<KdLeaf *>::const_iterator kit, kit_end = tObject->mKdLeaves.end();
1161                       
1162                        for (kit = tObject->mKdLeaves.begin(); kit != kit_end; ++ kit)
1163                        {
1164                                KdLeaf *node = dynamic_cast<KdLeaf *>(*kit);
1165
1166                                if (ci.type == SortableEntry::ERayMin)
1167                                {
1168                                        AddContriToPvs(node, pvsLeft);
1169                                }
1170                                else if (ci.type == SortableEntry::ERayMax)
1171                                {
1172                                        RemoveContriFromPvs(node, pvsRight);
1173                                }
1174                        }
1175                }
1176        }
1177}
1178
1179
1180float VspTree::EvalLocalCostHeuristics(const RayInfoContainer &rays,
1181                                                                           const AxisAlignedBox3 &box,
1182                                                                           int pvsSize,
1183                                                                           const int axis,
1184                                                                           float &position)
1185{
1186        const float minBox = box.Min(axis);
1187        const float maxBox = box.Max(axis);
1188
1189        const float sizeBox = maxBox - minBox;
1190
1191        const float minBand = minBox + mMinBand * sizeBox;
1192        const float maxBand = minBox + mMaxBand * sizeBox;
1193
1194        SortSplitCandidates(rays, axis, minBand, maxBand);
1195
1196        // prepare the sweep
1197        // (note: returns pvs size, so there would be no need
1198        // to give pvs size as argument)
1199        pvsSize = PrepareHeuristics(rays);
1200
1201        Debug << "here45 pvs: " << pvsSize << endl;
1202
1203        // go through the lists, count the number of objects left and right
1204        // and evaluate the following cost funcion:
1205        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
1206
1207        int pvsl = 0;
1208        int pvsr = pvsSize;
1209
1210        int pvsBack = pvsl;
1211        int pvsFront = pvsr;
1212
1213        float sum = (float)pvsSize * sizeBox;
1214        float minSum = 1e20f;
1215
1216       
1217        // if no good split can be found, take mid split
1218        position = minBox + 0.5f * sizeBox;
1219       
1220        // the relative cost ratio
1221        float ratio = 99999999.0f;
1222        bool splitPlaneFound = false;
1223
1224        Intersectable::NewMail();
1225        KdLeaf::NewMail();
1226
1227
1228        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
1229
1230Debug << "****************" << endl;
1231        //-- traverse through visibility events
1232
1233        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
1234        {
1235                EvalPvsIncr(*ci, pvsl, pvsr);
1236
1237                // Note: sufficient to compare size of bounding boxes of front and back side?
1238                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
1239                {
1240                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
1241
1242                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << "\t cost= " << sum << endl;
1243
1244                        if (sum < minSum)
1245                        {
1246                                splitPlaneFound = true;
1247
1248                                minSum = sum;
1249                                position = (*ci).value;
1250                               
1251                                pvsBack = pvsl;
1252                                pvsFront = pvsr;
1253                        }
1254                }
1255        }
1256       
1257       
1258        // -- compute cost
1259
1260        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1261        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1262
1263        const float pOverall = sizeBox;
1264        const float pBack = position - minBox;
1265        const float pFront = maxBox - position;
1266
1267        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
1268    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1269        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1270       
1271        const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
1272        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1273
1274        if (splitPlaneFound)
1275        {
1276                ratio = newRenderCost / oldRenderCost;
1277        }
1278       
1279        //if (axis != 1)
1280        Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1281              <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
1282
1283        return ratio;
1284}
1285
1286
1287float VspTree::SelectSplitPlane(const VspTraversalData &tData,
1288                                                                AxisAlignedPlane &plane,
1289                                                                float &pFront,
1290                                                                float &pBack)
1291{
1292        float nPosition[3];
1293        float nCostRatio[3];
1294        float nProbFront[3];
1295        float nProbBack[3];
1296
1297        // create bounding box of node geometry
1298        AxisAlignedBox3 box = tData.mBoundingBox;
1299               
1300        int sAxis = 0;
1301        int bestAxis = -1;
1302
1303        // if we use some kind of specialised fixed axis
1304    const bool useSpecialAxis =
1305                mOnlyDrivingAxis || mCirculatingAxis;
1306        //Debug << "data: " << tData.mBoundingBox << " pvs " << tData.mPvs << endl;
1307        if (mCirculatingAxis)
1308        {
1309                int parentAxis = 0;
1310                VspNode *parent = tData.mNode->GetParent();
1311
1312                if (parent)
1313                        parentAxis = dynamic_cast<VspInterior *>(parent)->GetAxis();
1314
1315                sAxis = (parentAxis + 1) % 3;
1316        }
1317        else if (mOnlyDrivingAxis)
1318        {
1319                sAxis = box.Size().DrivingAxis();
1320        }
1321        //sAxis = 2;
1322        for (int axis = 0; axis < 3; ++ axis)
1323        {
1324                if (!useSpecialAxis || (axis == sAxis))
1325                {
1326                        //-- place split plane using heuristics
1327
1328                        if (mUseCostHeuristics)
1329                        {
1330                                nCostRatio[axis] =
1331                                        EvalLocalCostHeuristics(*tData.mRays,
1332                                                                                        box,
1333                                                                                        tData.mPvs,
1334                                                                                        axis,
1335                                                                                        nPosition[axis]);                       
1336                        }
1337                        else //-- split plane position is spatial median
1338                        {
1339                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1340
1341                                nCostRatio[axis] = EvalLocalSplitCost(tData,
1342                                                                                                          box,
1343                                                                                                          axis,
1344                                                                                                          nPosition[axis],
1345                                                                                                          nProbFront[axis],
1346                                                                                                          nProbBack[axis]);
1347                        }
1348                                               
1349                        if (bestAxis == -1)
1350                        {
1351                                bestAxis = axis;
1352                        }
1353                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1354                        {
1355                                bestAxis = axis;
1356                        }
1357                }
1358        }
1359
1360
1361        //-- assign values
1362       
1363        plane.mAxis = bestAxis;
1364        // split plane position
1365    plane.mPosition = nPosition[bestAxis];
1366
1367        pFront = nProbFront[bestAxis];
1368        pBack = nProbBack[bestAxis];
1369
1370        //Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
1371        return nCostRatio[bestAxis];
1372}
1373
1374
1375float VspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
1376                                                                          const VspTraversalData &data) const
1377{
1378#if 0
1379        return (float)-data.mDepth;
1380#endif
1381        float pvsFront = 0;
1382        float pvsBack = 0;
1383        float totalPvs = 0;
1384
1385        // probability that view point lies in back / front node
1386        float pOverall = data.mProbability;
1387        float pFront = 0;
1388        float pBack = 0;
1389
1390
1391        // create unique ids for pvs heuristics
1392        Intersectable::NewMail();
1393       
1394        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1395
1396        for (rit = data.mRays->begin(); rit != rit_end; ++ rit)
1397        {
1398                RayInfo rayInf = *rit;
1399
1400                float t;
1401                VssRay *ray = rayInf.mRay;
1402
1403                const int cf =
1404                        rayInf.ComputeRayIntersection(candidatePlane.mAxis,
1405                                                                                  candidatePlane.mPosition, t);
1406
1407                // find front and back pvs for origing and termination object
1408                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1409                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1410        }
1411
1412
1413        AxisAlignedBox3 frontBox;
1414        AxisAlignedBox3 backBox;
1415
1416        data.mBoundingBox.Split(candidatePlane.mAxis, candidatePlane.mPosition, frontBox, backBox);
1417
1418        pFront = frontBox.GetVolume();
1419        pBack = pOverall - pFront;
1420               
1421
1422        //-- pvs rendering heuristics
1423        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1424        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1425
1426        //-- only render cost heuristics or combined with standard deviation
1427        const float penaltyOld = EvalPvsPenalty((int)totalPvs, lowerPvsLimit, upperPvsLimit);
1428    const float penaltyFront = EvalPvsPenalty((int)pvsFront, lowerPvsLimit, upperPvsLimit);
1429        const float penaltyBack = EvalPvsPenalty((int)pvsBack, lowerPvsLimit, upperPvsLimit);
1430                       
1431        const float oldRenderCost = pOverall * penaltyOld;
1432        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1433
1434        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
1435        const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBoundingBox.GetVolume();
1436       
1437        // take render cost of node into account
1438        // otherwise danger of being stuck in a local minimum!!
1439        const float factor = 0.99f;
1440
1441        const float normalizedOldRenderCost = oldRenderCost / mBoundingBox.GetVolume();
1442        return factor * renderCostDecrease + (1.0f - factor) * normalizedOldRenderCost;
1443}
1444
1445
1446float VspTree::EvalLocalSplitCost(const VspTraversalData &data,
1447                                                                  const AxisAlignedBox3 &box,
1448                                                                  const int axis,
1449                                                                  const float &position,
1450                                                                  float &pFront,
1451                                                                  float &pBack) const
1452{
1453        float pvsTotal = 0;
1454        float pvsFront = 0;
1455        float pvsBack = 0;
1456       
1457        // create unique ids for pvs heuristics
1458        Intersectable::NewMail();
1459
1460        const int pvsSize = data.mPvs;
1461
1462        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1463
1464        // this is the main ray classification loop!
1465        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
1466        {
1467                // determine the side of this ray with respect to the plane
1468                float t;
1469                const int side = (*rit).ComputeRayIntersection(axis, position, t);
1470       
1471                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
1472                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
1473        }
1474
1475        //-- pvs heuristics
1476        float pOverall;
1477
1478        //-- compute heurstics
1479       
1480        pOverall = data.mProbability;
1481        // we take simplified computation for mid split
1482        pBack = pFront = pOverall * 0.5f;
1483       
1484       
1485#ifdef _DEBUG
1486        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
1487        Debug << pFront << " " << pBack << " " << pOverall << endl;
1488#endif
1489
1490       
1491        const float newCost = pvsBack * pBack + pvsFront * pFront;
1492        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
1493
1494        return  (mCtDivCi + newCost) / oldCost;
1495}
1496
1497
1498void VspTree::AddObjToPvs(Intersectable *obj,
1499                                                  const int cf,
1500                                                  float &frontPvs,
1501                                                  float &backPvs,
1502                                                  float &totalPvs) const
1503{
1504        if (!obj)
1505                return;
1506
1507        //const float renderCost = mViewCellsManager->EvalRenderCost(obj);
1508        const int renderCost = 1;
1509
1510        // object in no pvs => new
1511        if (!obj->Mailed() && !obj->Mailed(1) && !obj->Mailed(2))
1512        {
1513                totalPvs += renderCost;
1514        }
1515
1516        // TODO: does this really belong to no pvs?
1517        //if (cf == Ray::COINCIDENT) return;
1518
1519        if (cf >= 0) // front pvs
1520        {
1521                if (!obj->Mailed() && !obj->Mailed(2))
1522                {
1523                        frontPvs += renderCost;
1524               
1525                        // already in back pvs => in both pvss
1526                        if (obj->Mailed(1))
1527                                obj->Mail(2);
1528                        else
1529                                obj->Mail();
1530                }
1531        }
1532
1533        if (cf <= 0) // back pvs
1534        {
1535                if (!obj->Mailed(1) && !obj->Mailed(2))
1536                {
1537                        backPvs += renderCost;
1538               
1539                        // already in front pvs => in both pvss
1540                        if (obj->Mailed())
1541                                obj->Mail(2);
1542                        else
1543                                obj->Mail(1);
1544                }
1545        }
1546}
1547
1548
1549void VspTree::CollectLeaves(vector<VspLeaf *> &leaves,
1550                                                           const bool onlyUnmailed,
1551                                                           const int maxPvsSize) const
1552{
1553        stack<VspNode *> nodeStack;
1554        nodeStack.push(mRoot);
1555
1556        while (!nodeStack.empty())
1557        {
1558                VspNode *node = nodeStack.top();
1559                nodeStack.pop();
1560               
1561                if (node->IsLeaf())
1562                {
1563                        // test if this leaf is in valid view space
1564                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1565                        if (leaf->TreeValid() &&
1566                                (!onlyUnmailed || !leaf->Mailed()) &&
1567                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
1568                        {
1569                                leaves.push_back(leaf);
1570                        }
1571                }
1572                else
1573                {
1574                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1575
1576                        nodeStack.push(interior->GetBack());
1577                        nodeStack.push(interior->GetFront());
1578                }
1579        }
1580}
1581
1582
1583AxisAlignedBox3 VspTree::GetBoundingBox() const
1584{
1585        return mBoundingBox;
1586}
1587
1588
1589VspNode *VspTree::GetRoot() const
1590{
1591        return mRoot;
1592}
1593
1594
1595void VspTree::EvaluateLeafStats(const VspTraversalData &data)
1596{
1597        // the node became a leaf -> evaluate stats for leafs
1598        VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);
1599
1600
1601        if (data.mPvs > mVspStats.maxPvs)
1602        {
1603                mVspStats.maxPvs = data.mPvs;
1604        }
1605
1606        mVspStats.pvs += data.mPvs;
1607
1608        if (data.mDepth < mVspStats.minDepth)
1609        {
1610                mVspStats.minDepth = data.mDepth;
1611        }
1612       
1613        if (data.mDepth >= mTermMaxDepth)
1614        {
1615        ++ mVspStats.maxDepthNodes;
1616                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
1617        }
1618
1619        // accumulate rays to compute rays /  leaf
1620        mVspStats.accumRays += (int)data.mRays->size();
1621
1622        if (data.mPvs < mTermMinPvs)
1623                ++ mVspStats.minPvsNodes;
1624
1625        if ((int)data.mRays->size() < mTermMinRays)
1626                ++ mVspStats.minRaysNodes;
1627
1628        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1629                ++ mVspStats.maxRayContribNodes;
1630
1631        if (data.mProbability <= mTermMinProbability)
1632                ++ mVspStats.minProbabilityNodes;
1633       
1634        // accumulate depth to compute average depth
1635        mVspStats.accumDepth += data.mDepth;
1636
1637        ++ mCreatedViewCells;
1638
1639#ifdef _DEBUG
1640        Debug << "BSP stats: "
1641                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1642                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1643                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1644                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "), "
1645                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1646#endif
1647}
1648
1649
1650void VspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
1651{
1652        ViewCell::NewMail();
1653        CollectViewCells(mRoot, onlyValid, viewCells, true);
1654}
1655
1656
1657void VspTree::CollapseViewCells()
1658{
1659// TODO
1660#if HAS_TO_BE_REDONE
1661        stack<VspNode *> nodeStack;
1662
1663        if (!mRoot)
1664                return;
1665
1666        nodeStack.push(mRoot);
1667       
1668        while (!nodeStack.empty())
1669        {
1670                VspNode *node = nodeStack.top();
1671                nodeStack.pop();
1672               
1673                if (node->IsLeaf())
1674        {
1675                        BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1676
1677                        if (!viewCell->GetValid())
1678                        {
1679                                BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1680       
1681                                ViewCellContainer leaves;
1682                                mViewCellsTree->CollectLeaves(viewCell, leaves);
1683
1684                                ViewCellContainer::const_iterator it, it_end = leaves.end();
1685
1686                                for (it = leaves.begin(); it != it_end; ++ it)
1687                                {
1688                                        VspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
1689                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
1690                                        ++ mVspStats.invalidLeaves;
1691                                }
1692
1693                                // add to unbounded view cell
1694                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
1695                                DEL_PTR(viewCell);
1696                        }
1697                }
1698                else
1699                {
1700                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1701               
1702                        nodeStack.push(interior->GetFront());
1703                        nodeStack.push(interior->GetBack());
1704                }
1705        }
1706
1707        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1708#endif
1709}
1710
1711
1712void VspTree::CollectRays(VssRayContainer &rays)
1713{
1714        vector<VspLeaf *> leaves;
1715
1716        vector<VspLeaf *>::const_iterator lit, lit_end = leaves.end();
1717
1718        for (lit = leaves.begin(); lit != lit_end; ++ lit)
1719        {
1720                VspLeaf *leaf = *lit;
1721                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
1722
1723                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
1724                        rays.push_back(*rit);
1725        }
1726}
1727
1728
1729void VspTree::SetViewCellsManager(ViewCellsManager *vcm)
1730{
1731        mViewCellsManager = vcm;
1732}
1733
1734
1735void VspTree::ValidateTree()
1736{
1737        mVspStats.invalidLeaves = 0;
1738
1739        stack<VspNode *> nodeStack;
1740
1741        if (!mRoot)
1742                return;
1743
1744        nodeStack.push(mRoot);
1745
1746        while (!nodeStack.empty())
1747        {
1748                VspNode *node = nodeStack.top();
1749                nodeStack.pop();
1750               
1751                if (node->IsLeaf())
1752                {
1753                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1754
1755                        if (!leaf->GetViewCell()->GetValid())
1756                                ++ mVspStats.invalidLeaves;
1757
1758                        // validity flags don't match => repair
1759                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
1760                        {
1761                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
1762                                PropagateUpValidity(leaf);
1763                        }
1764                }
1765                else
1766                {
1767                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1768               
1769                        nodeStack.push(interior->GetFront());
1770                        nodeStack.push(interior->GetBack());
1771                }
1772        }
1773
1774        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1775}
1776
1777
1778
1779void VspTree::CollectViewCells(VspNode *root,
1780                                                                  bool onlyValid,
1781                                                                  ViewCellContainer &viewCells,
1782                                                                  bool onlyUnmailed) const
1783{
1784        stack<VspNode *> nodeStack;
1785
1786        if (!root)
1787                return;
1788
1789        nodeStack.push(root);
1790       
1791        while (!nodeStack.empty())
1792        {
1793                VspNode *node = nodeStack.top();
1794                nodeStack.pop();
1795               
1796                if (node->IsLeaf())
1797                {
1798                        if (!onlyValid || node->TreeValid())
1799                        {
1800                                ViewCellLeaf *leafVc = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1801
1802                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
1803                                               
1804                                if (!onlyUnmailed || !viewCell->Mailed())
1805                                {
1806                                        viewCell->Mail();
1807                                        viewCells.push_back(viewCell);
1808                                }
1809                        }
1810                }
1811                else
1812                {
1813                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1814               
1815                        nodeStack.push(interior->GetFront());
1816                        nodeStack.push(interior->GetBack());
1817                }
1818        }
1819}
1820
1821
1822int VspTree::FindNeighbors(VspLeaf *n,
1823                                                   vector<VspLeaf *> &neighbors,
1824                                                   const bool onlyUnmailed) const
1825{
1826        stack<VspNode *> nodeStack;
1827        nodeStack.push(mRoot);
1828
1829        const AxisAlignedBox3 box = GetBBox(n);
1830
1831        while (!nodeStack.empty())
1832        {
1833                VspNode *node = nodeStack.top();
1834                nodeStack.pop();
1835
1836                if (node->IsLeaf())
1837                {
1838                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1839
1840                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
1841                                neighbors.push_back(leaf);
1842                }
1843                else
1844                {
1845                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1846                       
1847                        if (interior->GetPosition() > box.Max(interior->GetAxis()))
1848                                nodeStack.push(interior->GetBack());
1849                        else
1850                        {
1851                                if (interior->GetPosition() < box.Min(interior->GetAxis()))
1852                                        nodeStack.push(interior->GetFront());
1853                                else
1854                                {
1855                                        // random decision
1856                                        nodeStack.push(interior->GetBack());
1857                                        nodeStack.push(interior->GetFront());
1858                                }
1859                        }
1860                }
1861        }
1862
1863        return (int)neighbors.size();
1864}
1865
1866
1867// Find random neighbor which was not mailed
1868VspLeaf *VspTree::GetRandomLeaf(const Plane3 &plane)
1869{
1870        stack<VspNode *> nodeStack;
1871        nodeStack.push(mRoot);
1872 
1873        int mask = rand();
1874 
1875        while (!nodeStack.empty())
1876        {
1877                VspNode *node = nodeStack.top();
1878               
1879                nodeStack.pop();
1880               
1881                if (node->IsLeaf())
1882                {
1883                        return dynamic_cast<VspLeaf *>(node);
1884                }
1885                else
1886                {
1887                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1888                        VspNode *next;
1889                       
1890                        if (GetBBox(interior->GetBack()).Side(plane) < 0)
1891                        {
1892                                next = interior->GetFront();
1893                        }
1894            else
1895                        {
1896                                if (GetBBox(interior->GetFront()).Side(plane) < 0)
1897                                {
1898                                        next = interior->GetBack();
1899                                }
1900                                else
1901                                {
1902                                        // random decision
1903                                        if (mask & 1)
1904                                                next = interior->GetBack();
1905                                        else
1906                                                next = interior->GetFront();
1907                                        mask = mask >> 1;
1908                                }
1909                        }
1910                       
1911                        nodeStack.push(next);
1912                }
1913        }
1914 
1915        return NULL;
1916}
1917
1918
1919VspLeaf *VspTree::GetRandomLeaf(const bool onlyUnmailed)
1920{
1921        stack<VspNode *> nodeStack;
1922
1923        nodeStack.push(mRoot);
1924
1925        int mask = rand();
1926
1927        while (!nodeStack.empty())
1928        {
1929                VspNode *node = nodeStack.top();
1930                nodeStack.pop();
1931
1932                if (node->IsLeaf())
1933                {
1934                        if ( (!onlyUnmailed || !node->Mailed()) )
1935                                return dynamic_cast<VspLeaf *>(node);
1936                }
1937                else
1938                {
1939                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1940
1941                        // random decision
1942                        if (mask & 1)
1943                                nodeStack.push(interior->GetBack());
1944                        else
1945                                nodeStack.push(interior->GetFront());
1946
1947                        mask = mask >> 1;
1948                }
1949        }
1950
1951        return NULL;
1952}
1953
1954
1955void VspTree::CollectPvs(const RayInfoContainer &rays,
1956                                                 ObjectContainer &objects) const
1957{
1958        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1959
1960        Intersectable::NewMail();
1961
1962        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1963        {
1964                VssRay *ray = (*rit).mRay;
1965
1966                Intersectable *object;
1967                object = ray->mOriginObject;
1968
1969        if (object)
1970                {
1971                        if (!object->Mailed())
1972                        {
1973                                object->Mail();
1974                                objects.push_back(object);
1975                        }
1976                }
1977
1978                object = ray->mTerminationObject;
1979
1980                if (object)
1981                {
1982                        if (!object->Mailed())
1983                        {
1984                                object->Mail();
1985                                objects.push_back(object);
1986                        }
1987                }
1988        }
1989}
1990
1991
1992int VspTree::ComputePvsSize(const RayInfoContainer &rays) const
1993{
1994        int pvsSize = 0;
1995
1996        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1997
1998        Intersectable::NewMail();
1999
2000        for (rit = rays.begin(); rit != rays.end(); ++ rit)
2001        {
2002                VssRay *ray = (*rit).mRay;
2003
2004                if (ray->mOriginObject)
2005                {
2006                        if (!ray->mOriginObject->Mailed())
2007                        {
2008                                ray->mOriginObject->Mail();
2009                                ++ pvsSize;
2010                        }
2011                }
2012                if (ray->mTerminationObject)
2013                {
2014                        if (!ray->mTerminationObject->Mailed())
2015                        {
2016                                ray->mTerminationObject->Mail();
2017                                ++ pvsSize;
2018                        }
2019                }
2020        }
2021
2022        return pvsSize;
2023}
2024
2025
2026float VspTree::GetEpsilon() const
2027{
2028        return mEpsilon;
2029}
2030
2031
2032int VspTree::CastLineSegment(const Vector3 &origin,
2033                                                         const Vector3 &termination,
2034                             ViewCellContainer &viewcells)
2035{
2036        int hits = 0;
2037
2038        float mint = 0.0f, maxt = 1.0f;
2039        const Vector3 dir = termination - origin;
2040
2041        stack<LineTraversalData> tStack;
2042
2043        Intersectable::NewMail();
2044        ViewCell::NewMail();
2045
2046        Vector3 entp = origin;
2047        Vector3 extp = termination;
2048
2049        VspNode *node = mRoot;
2050        VspNode *farChild;
2051
2052        float position;
2053        int axis;
2054
2055        while (1)
2056        {
2057                if (!node->IsLeaf())
2058                {
2059                        VspInterior *in = dynamic_cast<VspInterior *>(node);
2060                        position = in->GetPosition();
2061                        axis = in->GetAxis();
2062
2063                        if (entp[axis] <= position)
2064                        {
2065                                if (extp[axis] <= position)
2066                                {
2067                                        node = in->GetBack();
2068                                        // cases N1,N2,N3,P5,Z2,Z3
2069                                        continue;
2070                                } else
2071                                {
2072                                        // case N4
2073                                        node = in->GetBack();
2074                                        farChild = in->GetFront();
2075                                }
2076                        }
2077                        else
2078                        {
2079                                if (position <= extp[axis])
2080                                {
2081                                        node = in->GetFront();
2082                                        // cases P1,P2,P3,N5,Z1
2083                                        continue;
2084                                }
2085                                else
2086                                {
2087                                        node = in->GetFront();
2088                                        farChild = in->GetBack();
2089                                        // case P4
2090                                }
2091                        }
2092
2093                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2094                        // case N4 or P4
2095                        const float tdist = (position - origin[axis]) / dir[axis];
2096                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2097
2098                        extp = origin + dir * tdist;
2099                        maxt = tdist;
2100                }
2101                else
2102                {
2103                        // compute intersection with all objects in this leaf
2104                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2105                        ViewCell *vc = leaf->GetViewCell();
2106
2107                        if (!vc->Mailed())
2108                        {
2109                                vc->Mail();
2110                                viewcells.push_back(vc);
2111                                ++ hits;
2112                        }
2113#if 0
2114                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2115#endif
2116                        // get the next node from the stack
2117                        if (tStack.empty())
2118                                break;
2119
2120                        entp = extp;
2121                        mint = maxt;
2122                       
2123                        LineTraversalData &s  = tStack.top();
2124                        node = s.mNode;
2125                        extp = s.mExitPoint;
2126                        maxt = s.mMaxT;
2127                        tStack.pop();
2128                }
2129        }
2130
2131        return hits;
2132}
2133
2134
2135int VspTree::CastRay(Ray &ray)
2136{
2137        int hits = 0;
2138
2139        stack<LineTraversalData> tStack;
2140        const Vector3 dir = ray.GetDir();
2141
2142        float maxt, mint;
2143
2144        if (!mBoundingBox.GetRaySegment(ray, mint, maxt))
2145                return 0;
2146
2147        Intersectable::NewMail();
2148        ViewCell::NewMail();
2149
2150        Vector3 entp = ray.Extrap(mint);
2151        Vector3 extp = ray.Extrap(maxt);
2152
2153        const Vector3 origin = entp;
2154
2155        VspNode *node = mRoot;
2156        VspNode *farChild = NULL;
2157
2158        float position;
2159        int axis;
2160
2161        while (1)
2162        {
2163                if (!node->IsLeaf())
2164                {
2165                        VspInterior *in = dynamic_cast<VspInterior *>(node);
2166                        position = in->GetPosition();
2167                        axis = in->GetAxis();
2168
2169                        if (entp[axis] <= position)
2170                        {
2171                                if (extp[axis] <= position)
2172                                {
2173                                        node = in->GetBack();
2174                                        // cases N1,N2,N3,P5,Z2,Z3
2175                                        continue;
2176                                }
2177                                else
2178                                {
2179                                        // case N4
2180                                        node = in->GetBack();
2181                                        farChild = in->GetFront();
2182                                }
2183                        }
2184                        else
2185                        {
2186                                if (position <= extp[axis])
2187                                {
2188                                        node = in->GetFront();
2189                                        // cases P1,P2,P3,N5,Z1
2190                                        continue;
2191                                }
2192                                else
2193                                {
2194                                        node = in->GetFront();
2195                                        farChild = in->GetBack();
2196                                        // case P4
2197                                }
2198                        }
2199
2200                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2201                        // case N4 or P4
2202                        const float tdist = (position - origin[axis]) / dir[axis];
2203                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2204                        extp = origin + dir * tdist;
2205                        maxt = tdist;
2206                }
2207                else
2208                {
2209                        // compute intersection with all objects in this leaf
2210                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2211                        ViewCell *vc = leaf->GetViewCell();
2212
2213                        if (!vc->Mailed())
2214                        {
2215                                vc->Mail();
2216                                // todo: add view cells to ray
2217                                ++ hits;
2218                        }
2219#if 0
2220                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2221#endif
2222                        // get the next node from the stack
2223                        if (tStack.empty())
2224                                break;
2225
2226                        entp = extp;
2227                        mint = maxt;
2228                       
2229                        LineTraversalData &s  = tStack.top();
2230                        node = s.mNode;
2231                        extp = s.mExitPoint;
2232                        maxt = s.mMaxT;
2233                        tStack.pop();
2234                }
2235        }
2236
2237        return hits;
2238}
2239
2240
2241ViewCell *VspTree::GetViewCell(const Vector3 &point, const bool active)
2242{
2243        if (mRoot == NULL)
2244                return NULL;
2245
2246        stack<VspNode *> nodeStack;
2247        nodeStack.push(mRoot);
2248 
2249        ViewCellLeaf *viewcell = NULL;
2250 
2251        while (!nodeStack.empty()) 
2252        {
2253                VspNode *node = nodeStack.top();
2254                nodeStack.pop();
2255       
2256                if (node->IsLeaf())
2257                {
2258                        viewcell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
2259                        break;
2260                }
2261                else   
2262                {       
2263                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2264     
2265                        // random decision
2266                        if (interior->GetPosition() - point[interior->GetAxis()] < 0)
2267                                nodeStack.push(interior->GetBack());
2268                        else
2269                                nodeStack.push(interior->GetFront());
2270                }
2271        }
2272 
2273        if (active)
2274                return mViewCellsTree->GetActiveViewCell(viewcell);
2275        else
2276                return viewcell;
2277}
2278
2279
2280bool VspTree::ViewPointValid(const Vector3 &viewPoint) const
2281{
2282        VspNode *node = mRoot;
2283
2284        while (1)
2285        {
2286                // early exit
2287                if (node->TreeValid())
2288                        return true;
2289
2290                if (node->IsLeaf())
2291                        return false;
2292                       
2293                VspInterior *in = dynamic_cast<VspInterior *>(node);
2294                                       
2295                if (in->GetPosition() - viewPoint[in->GetAxis()] <= 0)
2296                {
2297                        node = in->GetBack();
2298                }
2299                else
2300                {
2301                        node = in->GetFront();
2302                }
2303        }
2304
2305        // should never come here
2306        return false;
2307}
2308
2309
2310void VspTree::PropagateUpValidity(VspNode *node)
2311{
2312        const bool isValid = node->TreeValid();
2313
2314        // propagative up invalid flag until only invalid nodes exist over this node
2315        if (!isValid)
2316        {
2317                while (!node->IsRoot() && node->GetParent()->TreeValid())
2318                {
2319                        node = node->GetParent();
2320                        node->SetTreeValid(false);
2321                }
2322        }
2323        else
2324        {
2325                // propagative up valid flag until one of the subtrees is invalid
2326                while (!node->IsRoot() && !node->TreeValid())
2327                {
2328            node = node->GetParent();
2329                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2330                       
2331                        // the parent is valid iff both leaves are valid
2332                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
2333                                                           interior->GetFront()->TreeValid());
2334                }
2335        }
2336}
2337
2338#if ZIPPED_VIEWCELLS
2339bool VspTree::Export(ogzstream &stream)
2340#else
2341bool VspTree::Export(ofstream &stream)
2342#endif
2343{
2344        ExportNode(mRoot, stream);
2345
2346        return true;
2347}
2348
2349
2350#if ZIPPED_VIEWCELLS
2351void VspTree::ExportNode(VspNode *node, ogzstream &stream)
2352#else
2353void VspTree::ExportNode(VspNode *node, ofstream &stream)
2354#endif
2355{
2356        if (node->IsLeaf())
2357        {
2358                VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2359                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
2360
2361                int id = -1;
2362                if (viewCell != mOutOfBoundsCell)
2363                        id = viewCell->GetId();
2364
2365                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
2366        }
2367        else
2368        {       
2369                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2370       
2371                AxisAlignedPlane plane = interior->GetPlane();
2372                stream << "<Interior plane=\"" << plane.mPosition << " "
2373                           << plane.mAxis << "\">" << endl;
2374
2375                ExportNode(interior->GetBack(), stream);
2376                ExportNode(interior->GetFront(), stream);
2377
2378                stream << "</Interior>" << endl;
2379        }
2380}
2381
2382
2383int VspTree::SplitRays(const AxisAlignedPlane &plane,
2384                                           RayInfoContainer &rays,
2385                                           RayInfoContainer &frontRays,
2386                                           RayInfoContainer &backRays) const
2387{
2388        int splits = 0;
2389
2390        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2391
2392        for (rit = rays.begin(); rit != rit_end; ++ rit)
2393        {
2394                RayInfo bRay = *rit;
2395               
2396                VssRay *ray = bRay.mRay;
2397                float t;
2398
2399                // get classification and receive new t
2400                //-- test if start point behind or in front of plane
2401                const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
2402                       
2403
2404#if 1
2405                if (side == 0)
2406                {
2407                        ++ splits;
2408
2409                        if (ray->HasPosDir(plane.mAxis))
2410                        {
2411                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2412                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2413                        }
2414                        else
2415                        {
2416                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2417                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2418                        }
2419                }
2420                else if (side == 1)
2421                {
2422                        frontRays.push_back(bRay);
2423                }
2424                else
2425                {
2426                        backRays.push_back(bRay);
2427                }
2428#else
2429                if (side == 0)
2430                {
2431                        ++ splits;
2432
2433                        if (ray->HasPosDir(plane.mAxis))
2434                        {
2435                                backRays.push_back(RayInfo(ray, bRay.GetMaxT(), t));
2436                                frontRays.push_back(RayInfo(ray, t, bRay.GetMinT()));
2437                        }
2438                        else
2439                        {
2440                                frontRays.push_back(RayInfo(ray, bRay.GetMaxT(), t));
2441                                backRays.push_back(RayInfo(ray, t, bRay.GetMinT()));
2442                        }
2443                }
2444                else if (side == 1)
2445                {
2446                        backRays.push_back(bRay);
2447                }
2448                else
2449                {
2450                        frontRays.push_back(bRay);
2451                       
2452                }
2453#endif
2454        }
2455
2456        return splits;
2457}
2458
2459
2460AxisAlignedBox3 VspTree::GetBBox(VspNode *node) const
2461{
2462        if (!node->GetParent())
2463                return mBoundingBox;
2464
2465        if (!node->IsLeaf())
2466        {
2467                return (dynamic_cast<VspInterior *>(node))->GetBoundingBox();           
2468        }
2469
2470        VspInterior *parent = dynamic_cast<VspInterior *>(node->GetParent());
2471
2472        AxisAlignedBox3 box(parent->GetBoundingBox());
2473
2474        if (parent->GetFront() == node)
2475      box.SetMin(parent->GetAxis(), parent->GetPosition());
2476    else
2477      box.SetMax(parent->GetAxis(), parent->GetPosition());
2478
2479        return box;
2480}
2481
2482
2483
2484
2485int VspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box,
2486                                                                         ViewCellContainer &viewCells) const
2487{
2488        stack<VspNode *> nodeStack;
2489 
2490        ViewCell::NewMail();
2491
2492        while (!nodeStack.empty())
2493        {
2494                VspNode *node = nodeStack.top();
2495                nodeStack.pop();
2496
2497                const AxisAlignedBox3 bbox = GetBBox(node);
2498
2499                if (bbox.Includes(box))
2500                {
2501                        // node geometry is contained in box
2502                        CollectViewCells(node, true, viewCells, true);
2503                }
2504                else if (Overlap(bbox, box))
2505                {
2506                        if (node->IsLeaf())
2507                        {
2508                                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2509                       
2510                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
2511                                {
2512                                        leaf->GetViewCell()->Mail();
2513                                        viewCells.push_back(leaf->GetViewCell());
2514                                }
2515                        }
2516                        else
2517                        {
2518                                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2519                       
2520                                VspNode *first = interior->GetFront();
2521                                VspNode *second = interior->GetBack();
2522           
2523                                nodeStack.push(first);
2524                                nodeStack.push(second);
2525                        }
2526                }       
2527                // default: cull
2528        }
2529
2530        return (int)viewCells.size();
2531}
2532
2533
2534
2535/*****************************************************************/
2536/*                class OspTree implementation                   */
2537/*****************************************************************/
2538
2539
2540OspTree::OspTree():
2541mRoot(NULL)
2542#if TODO
2543mOutOfBoundsCell(NULL),
2544mStoreRays(false),
2545mUseRandomAxis(false),
2546mTimeStamp(1)
2547#endif
2548{
2549#if TODO
2550        bool randomize = false;
2551        Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
2552        if (randomize)
2553                Randomize(); // initialise random generator for heuristics
2554
2555        //-- termination criteria for autopartition
2556        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxDepth", mTermMaxDepth);
2557        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minPvs", mTermMinPvs);
2558        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minRays", mTermMinRays);
2559        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minProbability", mTermMinProbability);
2560        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxRayContribution", mTermMaxRayContribution);
2561       
2562        Environment::GetSingleton()->GetIntValue("VspTree.Termination.missTolerance", mTermMissTolerance);
2563        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxViewCells", mMaxViewCells);
2564
2565        //-- max cost ratio for early tree termination
2566        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxCostRatio", mTermMaxCostRatio);
2567
2568        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
2569        Environment::GetSingleton()->GetIntValue("VspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
2570
2571        // HACK//mTermMinPolygons = 25;
2572
2573        //-- factors for bsp tree split plane heuristics
2574        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.ct_div_ci", mCtDivCi);
2575
2576        //-- partition criteria
2577        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.epsilon", mEpsilon);
2578
2579        // if only the driving axis is used for axis aligned split
2580        Environment::GetSingleton()->GetBoolValue("VspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
2581       
2582        //Environment::GetSingleton()->GetFloatValue("VspTree.maxTotalMemory", mMaxTotalMemory);
2583        Environment::GetSingleton()->GetFloatValue("VspTree.maxStaticMemory", mMaxMemory);
2584
2585        Environment::GetSingleton()->GetBoolValue("VspTree.useCostHeuristics", mUseCostHeuristics);
2586        Environment::GetSingleton()->GetBoolValue("VspTree.simulateOctree", mCirculatingAxis);
2587
2588
2589        char subdivisionStatsLog[100];
2590        Environment::GetSingleton()->GetStringValue("VspTree.subdivisionStats", subdivisionStatsLog);
2591        mSubdivisionStats.open(subdivisionStatsLog);
2592
2593        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.minBand", mMinBand);
2594        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.maxBand", mMaxBand);
2595       
2596        mSplitBorder = 0.1f;
2597
2598
2599        //-- debug output
2600
2601        Debug << "******* OSP options ******** " << endl;
2602    Debug << "max depth: " << mTermMaxDepth << endl;
2603        Debug << "min PVS: " << mTermMinPvs << endl;
2604        Debug << "min probabiliy: " << mTermMinProbability << endl;
2605        Debug << "min rays: " << mTermMinRays << endl;
2606        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
2607        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
2608        Debug << "miss tolerance: " << mTermMissTolerance << endl;
2609        Debug << "max view cells: " << mMaxViewCells << endl;
2610        Debug << "randomize: " << randomize << endl;
2611
2612        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
2613        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
2614        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
2615        Debug << "max memory: " << mMaxMemory << endl;
2616        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
2617        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
2618       
2619        Debug << "circulating axis: " << mCirculatingAxis << endl;
2620        Debug << "minband: " << mMinBand << endl;
2621        Debug << "maxband: " << mMaxBand << endl;
2622       
2623
2624        mSplitCandidates = new vector<SortableEntry>;
2625
2626        Debug << endl;
2627#endif
2628}
2629
2630
2631
2632void OspTree::SplitObjects(const AxisAlignedPlane & splitPlane,
2633                                                   const ObjectContainer &objects,
2634                                                   ObjectContainer &back,
2635                                                   ObjectContainer &front)
2636{
2637        ObjectContainer::const_iterator oit, oit_end = objects.end();
2638
2639    for (oit = objects.begin(); oit != oit_end; ++ oit)
2640        {
2641                // determine the side of this ray with respect to the plane
2642                const AxisAlignedBox3 box = (*oit)->GetBox();
2643
2644                if (box.Max(splitPlane.mAxis) >= splitPlane.mPosition)
2645            front.push_back(*oit);
2646   
2647                if (box.Min(splitPlane.mAxis) < splitPlane.mPosition)
2648                        back.push_back(*oit);
2649#if TODO   
2650                mStat.objectRefs -= (int)objects.size();
2651                mStat.objectRefs += objectsBack + objectsFront;
2652#endif
2653        }
2654}
2655
2656
2657KdInterior *OspTree::SubdivideNode(KdLeaf *leaf,
2658                                                                   const AxisAlignedPlane &splitPlane,
2659                                                                   const AxisAlignedBox3 &box,
2660                                                                   AxisAlignedBox3 &backBBox,
2661                                                                   AxisAlignedBox3 &frontBBox)
2662{
2663#if TODO
2664    mSpatialStat.nodes += 2;
2665        mSpatialStat.splits[axis];
2666#endif
2667
2668        // add the new nodes to the tree
2669        KdInterior *node = new KdInterior(leaf->mParent);
2670
2671        const int axis = splitPlane.mAxis;
2672        const float position = splitPlane.mPosition;
2673
2674        node->mAxis = axis;
2675        node->mPosition = position;
2676        node->mBox = box;
2677
2678    backBBox = box;
2679        frontBBox = box;
2680 
2681        // first count ray sides
2682        int objectsBack = 0;
2683        int objectsFront = 0;
2684 
2685        backBBox.SetMax(axis, position);
2686        frontBBox.SetMin(axis, position);
2687
2688        ObjectContainer::const_iterator mi, mi_end = leaf->mObjects.end();
2689
2690        for ( mi = leaf->mObjects.begin(); mi != mi_end; ++ mi)
2691        {
2692                // determine the side of this ray with respect to the plane
2693                const AxisAlignedBox3 box = (*mi)->GetBox();
2694               
2695                if (box.Max(axis) > position)
2696                        ++ objectsFront;
2697   
2698                if (box.Min(axis) < position)
2699                        ++ objectsBack;
2700        }
2701
2702        KdLeaf *back = new KdLeaf(node, objectsBack);
2703        KdLeaf *front = new KdLeaf(node, objectsFront);
2704
2705        // replace a link from node's parent
2706        if (leaf->mParent)
2707                leaf->mParent->ReplaceChildLink(leaf, node);
2708
2709        // and setup child links
2710        node->SetupChildLinks(back, front);
2711
2712        SplitObjects(splitPlane, leaf->mObjects, back->mObjects, front->mObjects);
2713
2714        ProcessLeafObjects(back, leaf);
2715    ProcessLeafObjects(front, leaf);
2716 
2717        //delete leaf;
2718        return node;
2719}
2720
2721
2722KdNode *OspTree::Subdivide(SplitQueue &tQueue,
2723                                                   OspSplitCandidate &splitCandidate,
2724                                                   const bool globalCriteriaMet)
2725{
2726        OspTraversalData &tData = splitCandidate.mParentData;
2727
2728        KdNode *newNode = tData.mNode;
2729
2730        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
2731        {       
2732                OspTraversalData tFrontData;
2733                OspTraversalData tBackData;
2734
2735                //-- continue subdivision
2736               
2737                // create new interior node and two leaf node
2738                const AxisAlignedPlane splitPlane = splitCandidate.mSplitPlane;
2739               
2740                newNode = SubdivideNode(dynamic_cast<KdLeaf *>(newNode),
2741                                                                splitPlane,
2742                                                                tData.mBoundingBox,
2743                                                                tFrontData.mBoundingBox,
2744                                                                tBackData.mBoundingBox);
2745       
2746                const int maxCostMisses = splitCandidate.mMaxCostMisses;
2747
2748                // how often was max cost ratio missed in this branch?
2749                tFrontData.mMaxCostMisses = maxCostMisses;
2750                tBackData.mMaxCostMisses = maxCostMisses;
2751                       
2752                //-- push the new split candidates on the queue
2753                OspSplitCandidate *frontCandidate = new OspSplitCandidate();
2754                OspSplitCandidate *backCandidate = new OspSplitCandidate();
2755
2756                EvalSplitCandidate(tFrontData, *frontCandidate);
2757                EvalSplitCandidate(tBackData, *backCandidate);
2758
2759                tQueue.push(frontCandidate);
2760                tQueue.push(backCandidate);
2761
2762                // delete old leaf node
2763                DEL_PTR(tData.mNode);
2764        }
2765
2766
2767        //-- terminate traversal
2768        if (newNode->IsLeaf())
2769        {
2770                //KdLeaf *leaf = dynamic_cast<KdLeaf *>(newNode);
2771                EvaluateLeafStats(tData);               
2772        }
2773
2774        //-- cleanup
2775        tData.Clear();
2776       
2777        return newNode;
2778}
2779
2780
2781void OspTree::EvalSplitCandidate(OspTraversalData &tData,
2782                                                                 OspSplitCandidate &splitCandidate)
2783{
2784        float frontProb;
2785        float backProb;
2786       
2787        KdLeaf *leaf = dynamic_cast<KdLeaf *>(tData.mNode);
2788
2789        // compute locally best split plane
2790        const bool success =
2791                SelectSplitPlane(tData, splitCandidate.mSplitPlane,     frontProb, backProb);
2792
2793        //TODO
2794        // compute global decrease in render cost
2795        splitCandidate.mPriority = EvalRenderCostDecrease(splitCandidate.mSplitPlane, tData);
2796        splitCandidate.mParentData = tData;
2797        splitCandidate.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
2798}
2799
2800
2801bool OspTree::LocalTerminationCriteriaMet(const OspTraversalData &data) const
2802{
2803        // matt: TODO
2804        return true;
2805    /*          (((int)data.mRays->size() <= mTermMinRays) ||
2806                 (data.mPvs <= mTermMinPvs)   ||
2807                 (data.mProbability <= mTermMinProbability) ||
2808                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
2809                 (data.mDepth >= mTermMaxDepth));*/
2810}
2811
2812
2813bool OspTree::GlobalTerminationCriteriaMet(const OspTraversalData &data) const
2814{
2815        // matt: TODO
2816        return true;
2817                /*(mOutOfMemory ||
2818                (mVspStats.Leaves() >= mMaxViewCells) ||
2819                (mGlobalCostMisses >= mTermGlobalCostMissTolerance));*/
2820}
2821
2822
2823void OspTree::EvaluateLeafStats(const OspTraversalData &data)
2824{
2825#if TODO
2826        // the node became a leaf -> evaluate stats for leafs
2827        VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);
2828
2829        if (data.mPvs > mVspStats.maxPvs)
2830        {
2831                mVspStats.maxPvs = data.mPvs;
2832        }
2833
2834        mVspStats.pvs += data.mPvs;
2835
2836        if (data.mDepth < mVspStats.minDepth)
2837        {
2838                mVspStats.minDepth = data.mDepth;
2839        }
2840       
2841        if (data.mDepth >= mTermMaxDepth)
2842        {
2843        ++ mVspStats.maxDepthNodes;
2844                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
2845        }
2846
2847        // accumulate rays to compute rays /  leaf
2848        mVspStats.accumRays += (int)data.mRays->size();
2849
2850        if (data.mPvs < mTermMinPvs)
2851                ++ mVspStats.minPvsNodes;
2852
2853        if ((int)data.mRays->size() < mTermMinRays)
2854                ++ mVspStats.minRaysNodes;
2855
2856        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
2857                ++ mVspStats.maxRayContribNodes;
2858
2859        if (data.mProbability <= mTermMinProbability)
2860                ++ mVspStats.minProbabilityNodes;
2861       
2862        // accumulate depth to compute average depth
2863        mVspStats.accumDepth += data.mDepth;
2864
2865        ++ mCreatedViewCells;
2866
2867#ifdef _DEBUG
2868        Debug << "BSP stats: "
2869                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
2870                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
2871                  << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
2872                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
2873                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "), "
2874                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
2875#endif
2876
2877#endif
2878}
2879
2880
2881float OspTree::EvalLocalCostHeuristics(KdLeaf *node,
2882                                                                           const AxisAlignedBox3 &box,
2883                                                                           const int axis,
2884                                                                           float &position,
2885                                                                           int &objectsBack,
2886                                                                           int &objectsFront)
2887{
2888       
2889        SortSplitCandidates(node, axis);
2890 
2891        // go through the lists, count the number of objects left and right
2892        // and evaluate the following cost funcion:
2893        // C = ct_div_ci  + (ol + or)/queries
2894       
2895        int pvsSize = PrepareHeuristics(node->mObjects);;
2896        int pvsl = 0, pvsr = pvsSize;
2897 
2898        const float minBox = box.Min(axis);
2899        const float maxBox = box.Max(axis);
2900
2901        const float sizeBox = maxBox - minBox;
2902       
2903                       
2904        // if no good split can be found, take mid split
2905        position = minBox + 0.5f * sizeBox;
2906       
2907        // the relative cost ratio
2908        float ratio = 99999999.0f;
2909        bool splitPlaneFound = false;
2910 
2911        float minBand = minBox + mSplitBorder * (maxBox - minBox);
2912        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
2913
2914    float minSum = 1e20f;
2915
2916        int pvsBack = pvsl;
2917        int pvsFront = pvsr;
2918
2919        float sum = (float)pvsSize * sizeBox;
2920
2921        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
2922
2923        //-- traverse through visibility events
2924        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
2925        {
2926                EvalPvsIncr(*ci, pvsl, pvsr);
2927
2928                // Note: sufficient to compare size of bounding boxes of front and back side?
2929                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
2930                {
2931                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
2932
2933                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << "\t cost= " << sum << endl;
2934
2935                        if (sum < minSum)
2936                        {
2937                                splitPlaneFound = true;
2938
2939                                minSum = sum;
2940                                position = (*ci).value;
2941                               
2942                                pvsBack = pvsl;
2943                                pvsFront = pvsr;
2944                        }
2945                }
2946        }
2947       
2948       
2949        // -- compute cost
2950
2951        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
2952        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
2953
2954        const float pOverall = sizeBox;
2955        const float pBack = position - minBox;
2956        const float pFront = maxBox - position;
2957
2958        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
2959    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
2960        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
2961       
2962        const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
2963        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
2964
2965        if (splitPlaneFound)
2966        {
2967                ratio = newRenderCost / oldRenderCost;
2968        }
2969       
2970        //if (axis != 1)
2971        Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
2972              <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
2973
2974        return ratio;
2975}
2976
2977void OspTree::SortSplitCandidates(KdLeaf *node, const int axis)
2978{
2979        mSplitCandidates->clear();
2980
2981        int requestedSize = 2*(int)node->mObjects.size();
2982       
2983        // creates a sorted split candidates array
2984        if (mSplitCandidates->capacity() > 500000 &&
2985                requestedSize < (int)(mSplitCandidates->capacity()/10))
2986        {
2987                delete mSplitCandidates;
2988                mSplitCandidates = new vector<SortableEntry>;
2989        }
2990
2991        mSplitCandidates->reserve(requestedSize);
2992       
2993        ObjectContainer::const_iterator mi, mi_end = node->mObjects.end();
2994
2995        // insert all queries
2996        for(mi = node->mObjects.begin(); mi != mi_end; ++ mi)
2997        {
2998                AxisAlignedBox3 box = (*mi)->GetBox();
2999
3000                mSplitCandidates->push_back(SortableEntry(SortableEntry::BOX_MIN,
3001                                        box.Min(axis), *mi));
3002
3003
3004                mSplitCandidates->push_back(SortableEntry(SortableEntry::BOX_MAX,
3005                                        box.Max(axis), *mi));
3006        }
3007
3008        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
3009}
3010
3011
3012
3013int OspTree::PrepareHeuristics(Intersectable *object)
3014{
3015        ViewCellPvsMap::const_iterator vit, vit_end = object->mViewCellPvs.mEntries.end();
3016       
3017        int pvsSize = 0;
3018
3019        for (vit = object->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3020        {
3021        ViewCell *vc = (*vit).first;
3022
3023                if (!vc->Mailed())
3024                {
3025                        vc->Mail();
3026                        vc->mCounter = 1;
3027                        ++ pvsSize;
3028                }
3029                else
3030                {
3031                        ++ vc->mCounter;
3032                }
3033        }
3034
3035        return pvsSize;
3036}
3037
3038
3039int OspTree::PrepareHeuristics(const ObjectContainer &objects)
3040{       
3041        Intersectable::NewMail();
3042        ViewCell::NewMail();
3043
3044        int pvsSize = 0;
3045
3046        ObjectContainer::const_iterator oit, oit_end = objects.end();
3047
3048        //-- set all pvs as belonging to the front pvs
3049        for (oit = objects.begin(); oit != oit_end; ++ oit)
3050        {
3051                Intersectable *obj = *oit;
3052
3053                pvsSize += PrepareHeuristics(obj);             
3054        }
3055
3056        return pvsSize;
3057}
3058
3059
3060void OspTree::EvalPvsIncr(const SortableEntry &ci,
3061                                                  int &pvsLeft,
3062                                                  int &pvsRight) const
3063{
3064        Intersectable *obj = ci.mObject;
3065
3066        switch (ci.type)
3067        {
3068                case SortableEntry::BOX_MIN:
3069                        AddContriToPvs(obj, pvsLeft);
3070                        break;
3071                       
3072                case SortableEntry::BOX_MAX:
3073                        RemoveContriFromPvs(obj, pvsRight);
3074                        break;
3075        }
3076}
3077
3078
3079void OspTree::RemoveContriFromPvs(Intersectable *object, int &pvs) const
3080{
3081        ViewCellPvsMap::const_iterator vit, vit_end = object->mViewCellPvs.mEntries.end();
3082       
3083        for (vit = object->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3084        {
3085        ViewCell *vc = (*vit).first;
3086
3087                if (-- vc->mCounter == 0)
3088                {
3089                        -- pvs;
3090                }
3091        }
3092}
3093
3094
3095void OspTree::AddContriToPvs(Intersectable *object, int &pvs) const
3096{
3097        ViewCellPvsMap::const_iterator vit, vit_end = object->mViewCellPvs.mEntries.end();
3098
3099        for (vit = object->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3100        {
3101        ViewCell *vc = (*vit).first;
3102
3103                if (!vc->Mailed())
3104                {
3105                        vc->Mail();             
3106                        ++ pvs;
3107                }
3108        }
3109}
3110
3111
3112float OspTree::SelectSplitPlane(const OspTraversalData &tData,
3113                                                                AxisAlignedPlane &plane,
3114                                                                float &pFront,
3115                                                                float &pBack)
3116{
3117        float nPosition[3];
3118        float nCostRatio[3];
3119        float nProbFront[3];
3120        float nProbBack[3];
3121
3122        // create bounding box of node geometry
3123        AxisAlignedBox3 box = tData.mBoundingBox;
3124               
3125        int sAxis = 0;
3126        int bestAxis = -1;
3127
3128        if (mOnlyDrivingAxis)
3129        {
3130                sAxis = box.Size().DrivingAxis();
3131        }
3132
3133/*
3134        //sAxis = 2;
3135        for (int axis = 0; axis < 3; ++ axis)
3136        {
3137                if (!mOnlyDrivingAxis || (axis == sAxis))
3138                {
3139                        //-- place split plane using heuristics
3140
3141                        if (mUseCostHeuristics)
3142                        {
3143                                nCostRatio[axis] =
3144                                        EvalLocalCostHeuristics(*tData.mRays,
3145                                                                                        box,
3146                                                                                        tData.mPvs,
3147                                                                                        axis,
3148                                                                                        nPosition[axis]);                       
3149                        }
3150                        else //-- split plane position is spatial median
3151                        {
3152                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
3153
3154                                nCostRatio[axis] = EvalLocalSplitCost(tData,
3155                                                                                                          box,
3156                                                                                                          axis,
3157                                                                                                          nPosition[axis],
3158                                                                                                          nProbFront[axis],
3159                                                                                                          nProbBack[axis]);
3160                        }
3161                                               
3162                        if (bestAxis == -1)
3163                        {
3164                                bestAxis = axis;
3165                        }
3166                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
3167                        {
3168                                bestAxis = axis;
3169                        }
3170                }
3171        }
3172
3173*/
3174        //-- assign values
3175       
3176        plane.mAxis = bestAxis;
3177        // split plane position
3178    plane.mPosition = nPosition[bestAxis];
3179
3180        pFront = nProbFront[bestAxis];
3181        pBack = nProbBack[bestAxis];
3182
3183        //Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
3184        return nCostRatio[bestAxis];
3185}
3186
3187
3188float OspTree::EvalViewCellPvsIncr(Intersectable *object) const
3189{
3190        return 0;
3191}
3192
3193
3194float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
3195                                                                          const OspTraversalData &data) const
3196{
3197#if 0
3198        return (float)-data.mDepth;
3199#endif
3200
3201        float pvsFront = 0;
3202        float pvsBack = 0;
3203        float totalPvs = 0;
3204
3205        // probability that view point lies in back / front node
3206        float pOverall = data.mProbability;
3207        float pFront = 0;
3208        float pBack = 0;
3209
3210
3211        Intersectable::NewMail();
3212        ViewCell::NewMail();
3213
3214        KdLeaf *leaf = dynamic_cast<KdLeaf *>(data.mNode);
3215        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
3216
3217        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
3218        {
3219                Intersectable *obj = *oit;
3220                const AxisAlignedBox3 box = obj->GetBox();
3221
3222                if (box.Max(candidatePlane.mAxis) > candidatePlane.mPosition)
3223                        pvsFront += EvalViewCellPvsIncr(obj);
3224                if (box.Min(candidatePlane.mAxis) > candidatePlane.mPosition)
3225                        pvsBack += EvalViewCellPvsIncr(obj);
3226        }
3227
3228
3229        AxisAlignedBox3 frontBox;
3230        AxisAlignedBox3 backBox;
3231
3232        data.mBoundingBox.Split(candidatePlane.mAxis, candidatePlane.mPosition, frontBox, backBox);
3233
3234        pFront = frontBox.GetVolume();
3235        pBack = pOverall - pFront;
3236               
3237
3238        //-- pvs rendering heuristics
3239        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
3240        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
3241
3242        //-- only render cost heuristics or combined with standard deviation
3243        const float penaltyOld = EvalPvsPenalty((int)totalPvs, lowerPvsLimit, upperPvsLimit);
3244    const float penaltyFront = EvalPvsPenalty((int)pvsFront, lowerPvsLimit, upperPvsLimit);
3245        const float penaltyBack = EvalPvsPenalty((int)pvsBack, lowerPvsLimit, upperPvsLimit);
3246                       
3247        const float oldRenderCost = pOverall * penaltyOld;
3248        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
3249
3250        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
3251        const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBoundingBox.GetVolume();
3252       
3253        // take render cost of node into account
3254        // otherwise danger of being stuck in a local minimum!!
3255        const float factor = 0.99f;
3256
3257        const float normalizedOldRenderCost = oldRenderCost / mBoundingBox.GetVolume();
3258        return factor * renderCostDecrease + (1.0f - factor) * normalizedOldRenderCost;
3259}
3260
3261
3262void OspTree::PrepareConstruction(const ObjectContainer &objects,
3263                                                                  AxisAlignedBox3 *forcedBoundingBox)
3264{
3265        mOspStats.nodes = 1;
3266       
3267        if (forcedBoundingBox)
3268        {
3269                mBoundingBox = *forcedBoundingBox;
3270        }
3271        else // compute vsp tree bounding box
3272        {
3273                mBoundingBox.Initialize();
3274
3275                ObjectContainer::const_iterator oit, oit_end = objects.end();
3276
3277                //-- compute bounding box
3278        for (oit = objects.begin(); oit != oit_end; ++ oit)
3279                {
3280                        Intersectable *obj = *oit;
3281
3282                        // compute bounding box of view space
3283                        mBoundingBox.Include(obj->GetBox());
3284                        mBoundingBox.Include(obj->GetBox());
3285                }
3286
3287                mTermMinProbability *= mBoundingBox.GetVolume();
3288                mGlobalCostMisses = 0;
3289        }
3290}
3291
3292
3293void OspTree::ProcessLeafObjects(KdLeaf *leaf, KdLeaf *parent) const
3294{
3295        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
3296
3297        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
3298        {
3299                Intersectable *object = *oit;
3300
3301                if (parent)
3302                {
3303                        set<KdLeaf *>::iterator kdit = object->mKdLeaves.find(parent);
3304
3305                        // remove parent leaf
3306                        if (kdit != object->mKdLeaves.end())
3307                                object->mKdLeaves.erase(kdit);
3308                }
3309
3310                object->mKdLeaves.insert(leaf);
3311
3312                if (object->mKdLeaves.size() > 1)
3313                        leaf->mMultipleObjects.push_back(object);
3314        }
3315}
3316
3317
3318
3319/*********************************************************************/
3320/*                class HierarchyManager implementation              */
3321/*********************************************************************/
3322
3323
3324
3325HierarchyManager::HierarchyManager(VspTree &vspTree, OspTree &ospTree):
3326mVspTree(vspTree), mOspTree(ospTree)
3327{
3328}
3329
3330
3331SplitCandidate *HierarchyManager::NextSplitCandidate()
3332{
3333        SplitCandidate *splitCandidate = mTQueue.top();
3334        //Debug << "priority: " << splitCandidate->GetPriority() << endl;
3335        mTQueue.pop();
3336
3337        return splitCandidate;
3338}
3339
3340
3341void HierarchyManager::PrepareConstruction(const VssRayContainer &sampleRays,
3342                                                                                   const ObjectContainer &objects,
3343                                                                                   AxisAlignedBox3 *forcedViewSpace,
3344                                                                                   RayInfoContainer &rays)
3345{
3346        mVspTree.PrepareConstruction(sampleRays, forcedViewSpace);
3347
3348        long startTime = GetTime();
3349
3350        cout << "storing rays ... ";
3351
3352        Intersectable::NewMail();
3353
3354        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
3355
3356        //-- store rays
3357        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
3358        {
3359                VssRay *ray = *rit;
3360
3361                float minT, maxT;
3362
3363                static Ray hray;
3364                hray.Init(*ray);
3365               
3366                // TODO: not very efficient to implictly cast between rays types
3367                if (mVspTree.GetBoundingBox().GetRaySegment(hray, minT, maxT))
3368                {
3369                        float len = ray->Length();
3370
3371                        if (!len)
3372                                len = Limits::Small;
3373
3374                        rays.push_back(RayInfo(ray, minT / len, maxT / len));
3375                }
3376        }
3377
3378        cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
3379
3380
3381        int pvsSize = mVspTree.ComputePvsSize(rays);
3382
3383        // -- prepare view space partition
3384
3385        // add first candidate for view space partition
3386        mVspTree.mRoot = new VspLeaf();
3387        const float prop = mVspTree.mBoundingBox.GetVolume();
3388
3389        // first vsp traversal data
3390        VspTree::VspTraversalData vData(mVspTree.mRoot,
3391                                                                        0,
3392                                                                        &rays,
3393                                                                        //(int)objects.size(),
3394                                                                        pvsSize,       
3395                                                                        prop,
3396                                                                        mVspTree.mBoundingBox);
3397
3398
3399        // compute first split candidate
3400        VspTree::VspSplitCandidate *splitCandidate = new VspTree::VspSplitCandidate();
3401    mVspTree.EvalSplitCandidate(vData, *splitCandidate);
3402
3403        mTQueue.push(splitCandidate);
3404
3405
3406        //-- object space partition
3407
3408        mOspTree.PrepareConstruction(objects, forcedViewSpace);
3409
3410        // add first candidate for view space partition
3411        KdLeaf *leaf = new KdLeaf(NULL, 0);
3412        leaf->mObjects = objects;
3413
3414        mOspTree.mRoot = leaf;
3415       
3416
3417        // first osp traversal data
3418        OspTree::OspTraversalData oData(mOspTree.mRoot,
3419                                                                        0,
3420                                                                        &rays,
3421                                                                        pvsSize,
3422                                                                        prop,
3423                                                                        mOspTree.mBoundingBox);
3424
3425               
3426        mOspTree.ProcessLeafObjects(leaf, NULL);
3427
3428        // compute first split candidate
3429        OspTree::OspSplitCandidate *oSplitCandidate = new OspTree::OspSplitCandidate();
3430    mOspTree.EvalSplitCandidate(oData, *oSplitCandidate);
3431
3432        mTQueue.push(splitCandidate);
3433}
3434
3435
3436bool HierarchyManager::GlobalTerminationCriteriaMet(SplitCandidate *candidate) const
3437{
3438        if (candidate->Type() == SplitCandidate::VIEW_SPACE)
3439        {
3440                VspTree::VspSplitCandidate *sc =
3441                        dynamic_cast<VspTree::VspSplitCandidate *>(candidate);
3442               
3443                return mVspTree.GlobalTerminationCriteriaMet(sc->mParentData);
3444        }
3445        else
3446        {
3447                OspTree::OspSplitCandidate *sc =
3448                        dynamic_cast<OspTree::OspSplitCandidate *>(candidate);
3449               
3450                return mOspTree.GlobalTerminationCriteriaMet(sc->mParentData);
3451        }
3452
3453        return true;
3454}
3455
3456
3457void HierarchyManager::Construct(const VssRayContainer &sampleRays,
3458                                                                 const ObjectContainer &objects,
3459                                                                 AxisAlignedBox3 *forcedViewSpace)
3460{
3461        RayInfoContainer *rays = new RayInfoContainer();
3462       
3463        // prepare vsp and osp trees for traversal
3464        PrepareConstruction(sampleRays, objects, forcedViewSpace, *rays);
3465
3466        mVspTree.mVspStats.Reset();
3467        mVspTree.mVspStats.Start();
3468
3469        cout << "Constructing view space / object space tree ... \n";
3470        const long startTime = GetTime();       
3471       
3472        int i = 0;
3473        while (!FinishedConstruction())
3474        {
3475                SplitCandidate *splitCandidate = NextSplitCandidate();
3476           
3477                const bool globalTerminationCriteriaMet =
3478                        GlobalTerminationCriteriaMet(splitCandidate);
3479
3480                cout << "view cells: " << i ++ << endl;
3481
3482                // cost ratio of cost decrease / totalCost
3483                const float costRatio = splitCandidate->GetPriority() / mTotalCost;
3484                //Debug << "cost ratio: " << costRatio << endl;
3485
3486                if (costRatio < mTermMinGlobalCostRatio)
3487                        ++ mGlobalCostMisses;
3488
3489                //-- subdivide leaf node
3490                //-- either a object space or view space split
3491                if (splitCandidate->Type() == SplitCandidate::VIEW_SPACE)
3492                {
3493                        VspTree::VspSplitCandidate *sc =
3494                                dynamic_cast<VspTree::VspSplitCandidate *>(splitCandidate);
3495
3496                        VspNode *r = mVspTree.Subdivide(mTQueue, *sc, globalTerminationCriteriaMet);
3497                }
3498                else // object space split
3499                {                       
3500                        OspTree::OspSplitCandidate *sc =
3501                                dynamic_cast<OspTree::OspSplitCandidate *>(splitCandidate);
3502
3503                        KdNode *r = mOspTree.Subdivide(mTQueue, *sc, globalTerminationCriteriaMet);
3504                }
3505
3506                DEL_PTR(splitCandidate);
3507        }
3508
3509        cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl;
3510
3511        mVspTree.mVspStats.Stop();
3512}
3513
3514
3515bool HierarchyManager::FinishedConstruction()
3516{
3517        return mTQueue.empty();
3518}
3519
3520
3521void HierarchyManager::RepairQueue()
3522{
3523        // TODO
3524        // for each update of the view space partition:
3525        // the candidates from object space partition which
3526        // have been afected by the view space split (the kd split candidates
3527        // which saw the view cell which was split) must be reevaluated
3528        // (maybe not locally, just reinsert them into the queue)
3529        //
3530        // vice versa for the view cells
3531        // for each update of the object space partition
3532        // reevaluate split candidate for view cells which saw the split kd cell
3533        //
3534        // the priority queue update can be solved by implementing a binary heap
3535        // (explicit data structure, binary tree)
3536        // *) inserting and removal is efficient
3537        // *) search is not efficient => store pointer to queue element with each
3538        // split candidate
3539
3540        vector<SplitCandidate *> candidates;
3541
3542        while (!mTQueue.empty())
3543        {
3544                SplitCandidate *candidate = mTQueue.top();
3545                mTQueue.pop();
3546
3547                candidates.push_back(candidate);
3548        }
3549
3550        // Reinsert
3551
3552}
3553
3554
3555
3556/********************************************************/
3557/*             SplitHeap implementation                 */
3558/********************************************************/
3559
3560
3561SplitHeap::SplitHeap():mRoot(NULL)
3562{}
3563
3564void SplitHeap::Push(SplitCandidate *candidate)
3565{
3566        InsertTail(candidate);
3567
3568        // Swap until heap constaints fullfilled
3569        while (HeapViolated(candidate))
3570        {
3571                Swap(candidate, candidate->mParent);
3572        }
3573}
3574
3575
3576void SplitHeap::InsertTail(SplitCandidate *candidate)
3577{
3578}
3579
3580
3581bool SplitHeap::HeapViolated(SplitCandidate *candidate)
3582{
3583        return true;
3584}
3585       
3586SplitCandidate *SplitHeap::Pop()
3587{
3588
3589        return mRoot;
3590}
3591
3592void SplitHeap::Remove(SplitCandidate *candidate)
3593{
3594}
3595
3596}
Note: See TracBrowser for help on using the repository browser.