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

Revision 1149, 113.6 KB checked in by mattausch, 18 years ago (diff)
Line 
1#include <stack>
2#include <time.h>
3#include <iomanip>
4
5#include "ViewCell.h"
6#include "Plane3.h"
7#include "VspOspTree.h"
8#include "Mesh.h"
9#include "common.h"
10#include "Environment.h"
11#include "Polygon3.h"
12#include "Ray.h"
13#include "AxisAlignedBox3.h"
14#include "Exporter.h"
15#include "Plane3.h"
16#include "ViewCellsManager.h"
17#include "Beam.h"
18#include "KdTree.h"
19#include "KdIntersectable.h"
20
21
22namespace GtpVisibilityPreprocessor {
23
24#define USE_FIXEDPOINT_T 0
25
26
27//-- static members
28
29VspTree *VspTree::VspSplitCandidate::sVspTree = NULL;
30OspTree *OspTree::OspSplitCandidate::sOspTree = NULL;
31
32// pvs penalty can be different from pvs size
33inline static float EvalPvsPenalty(const int pvs,
34                                                                   const int lower,
35                                                                   const int upper)
36{
37        // clamp to minmax values
38        if (pvs < lower)
39                return (float)lower;
40        else if (pvs > upper)
41                return (float)upper;
42
43        return (float)pvs;
44}
45
46
47int VspNode::sMailId = 1;
48
49
50void VspTreeStatistics::Print(ostream &app) const
51{
52        app << "=========== VspTree statistics ===============\n";
53
54        app << setprecision(4);
55
56        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
57
58        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
59
60        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
61
62        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
63
64        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits[0] + splits[1] + splits[2] << endl;
65
66        app << "#N_SPLITS ( Number of splits in axes x y z)\n";
67
68        for (int i = 0; i < 3; ++ i)
69                app << splits[i] << " ";
70        app << endl;
71
72        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
73                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
74
75        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n"
76                << minPvsNodes * 100 / (double)Leaves() << endl;
77
78        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n"
79                << minRaysNodes * 100 / (double)Leaves() << endl;
80
81        app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
82                << maxCostNodes * 100 / (double)Leaves() << endl;
83
84        app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
85                << minProbabilityNodes * 100 / (double)Leaves() << endl;
86
87        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"
88                <<      maxRayContribNodes * 100 / (double)Leaves() << endl;
89
90        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
91
92        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
93
94        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
95
96        app << "#N_INVALIDLEAVES (number of invalid leaves )\n" << invalidLeaves << endl;
97
98        app << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
99        //app << "#N_PVS: " << pvs << endl;
100
101        app << "========== END OF VspTree statistics ==========\n";
102}
103
104
105void OspTreeStatistics::Print(ostream &app) const
106{
107        app << "=========== OspTree statistics ===============\n";
108
109        app << setprecision(4);
110
111        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
112
113        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
114
115        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
116
117        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
118
119        app << "#AXIS_ALIGNED_SPLITS (number of axis aligned splits)\n" << splits[0] + splits[1] + splits[2] << endl;
120
121        app << "#N_SPLITS ( Number of splits in axes x y z)\n";
122
123        for (int i = 0; i < 3; ++ i)
124                app << splits[i] << " ";
125        app << endl;
126
127        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
128                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
129
130        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n"
131                << minPvsNodes * 100 / (double)Leaves() << endl;
132
133        //app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n"
134        //      << minRaysNodes * 100 / (double)Leaves() << endl;
135
136        app << "#N_MAXCOSTNODES  ( Percentage of leaves with terminated because of max cost ratio )\n"
137                << maxCostNodes * 100 / (double)Leaves() << endl;
138
139        app << "#N_PMINPROBABILITYLEAVES  ( Percentage of leaves with mininum probability )\n"
140                << minProbabilityNodes * 100 / (double)Leaves() << endl;
141
142//      app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"
143//              <<      maxRayContribNodes * 100 / (double)Leaves() << endl;
144
145        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
146
147        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
148
149        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
150
151        app << "#N_INVALIDLEAVES (number of invalid leaves )\n" << invalidLeaves << endl;
152
153//      app << "#N_RAYS (number of rays / leaf)\n" << AvgRays() << endl;
154        //app << "#N_PVS: " << pvs << endl;
155
156        app << "========== END OF VspTree statistics ==========\n";
157}
158
159
160/******************************************************************/
161/*                  class VspNode implementation                  */
162/******************************************************************/
163
164
165VspNode::VspNode():
166mParent(NULL), mTreeValid(true), mTimeStamp(0)
167{}
168
169
170VspNode::VspNode(VspInterior *parent):
171mParent(parent), mTreeValid(true)
172{}
173
174
175bool VspNode::IsRoot() const
176{
177        return mParent == NULL;
178}
179
180
181VspInterior *VspNode::GetParent()
182{
183        return mParent;
184}
185
186
187void VspNode::SetParent(VspInterior *parent)
188{
189        mParent = parent;
190}
191
192
193bool VspNode::IsSibling(VspNode *n) const
194{
195        return  ((this != n) && mParent &&
196                         (mParent->GetFront() == n) || (mParent->GetBack() == n));
197}
198
199
200int VspNode::GetDepth() const
201{
202        int depth = 0;
203        VspNode *p = mParent;
204       
205        while (p)
206        {
207                p = p->mParent;
208                ++ depth;
209        }
210
211        return depth;
212}
213
214
215bool VspNode::TreeValid() const
216{
217        return mTreeValid;
218}
219
220
221void VspNode::SetTreeValid(const bool v)
222{
223        mTreeValid = v;
224}
225
226
227/****************************************************************/
228/*              class VspInterior implementation                */
229/****************************************************************/
230
231
232VspInterior::VspInterior(const AxisAlignedPlane &plane):
233mPlane(plane), mFront(NULL), mBack(NULL)
234{}
235
236
237VspInterior::~VspInterior()
238{
239        DEL_PTR(mFront);
240        DEL_PTR(mBack);
241}
242
243
244bool VspInterior::IsLeaf() const
245{
246        return false;
247}
248
249
250VspNode *VspInterior::GetBack()
251{
252        return mBack;
253}
254
255
256VspNode *VspInterior::GetFront()
257{
258        return mFront;
259}
260
261
262AxisAlignedPlane VspInterior::GetPlane() const
263{
264        return mPlane;
265}
266
267
268float VspInterior::GetPosition() const
269{
270        return mPlane.mPosition;
271}
272
273
274int VspInterior::GetAxis() const
275{
276        return mPlane.mAxis;
277}
278
279
280void VspInterior::ReplaceChildLink(VspNode *oldChild, VspNode *newChild)
281{
282        if (mBack == oldChild)
283                mBack = newChild;
284        else
285                mFront = newChild;
286}
287
288
289void VspInterior::SetupChildLinks(VspNode *front, VspNode *back)
290{
291    mBack = back;
292    mFront = front;
293}
294
295
296AxisAlignedBox3 VspInterior::GetBoundingBox() const
297{
298        return mBoundingBox;
299}
300
301
302void VspInterior::SetBoundingBox(const AxisAlignedBox3 &box)
303{
304        mBoundingBox = box;
305}
306
307
308int VspInterior::Type() const
309{
310        return Interior;
311}
312
313
314
315/****************************************************************/
316/*                  class VspLeaf implementation                */
317/****************************************************************/
318
319
320VspLeaf::VspLeaf(): mViewCell(NULL), mPvs(NULL)
321{
322}
323
324
325VspLeaf::~VspLeaf()
326{
327        DEL_PTR(mPvs);
328        CLEAR_CONTAINER(mVssRays);
329}
330
331
332int VspLeaf::Type() const
333{
334        return Leaf;
335}
336
337
338VspLeaf::VspLeaf(ViewCellLeaf *viewCell):
339mViewCell(viewCell)
340{
341}
342
343
344VspLeaf::VspLeaf(VspInterior *parent):
345VspNode(parent), mViewCell(NULL), mPvs(NULL)
346{}
347
348
349
350VspLeaf::VspLeaf(VspInterior *parent, ViewCellLeaf *viewCell):
351VspNode(parent), mViewCell(viewCell), mPvs(NULL)
352{
353}
354
355ViewCellLeaf *VspLeaf::GetViewCell() const
356{
357        return mViewCell;
358}
359
360void VspLeaf::SetViewCell(ViewCellLeaf *viewCell)
361{
362        mViewCell = viewCell;
363}
364
365
366bool VspLeaf::IsLeaf() const
367{
368        return true;
369}
370
371
372/******************************************************************************/
373/*                       class VspTree implementation                      */
374/******************************************************************************/
375
376
377VspTree::VspTree():
378mRoot(NULL),
379mOutOfBoundsCell(NULL),
380mStoreRays(false),
381mTimeStamp(1),
382mOspTree(NULL)
383{
384        bool randomize = false;
385        Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
386        if (randomize)
387                Randomize(); // initialise random generator for heuristics
388
389        //-- termination criteria for autopartition
390        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxDepth", mTermMaxDepth);
391        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minPvs", mTermMinPvs);
392        Environment::GetSingleton()->GetIntValue("VspTree.Termination.minRays", mTermMinRays);
393        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minProbability", mTermMinProbability);
394        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxRayContribution", mTermMaxRayContribution);
395       
396        Environment::GetSingleton()->GetIntValue("VspTree.Termination.missTolerance", mTermMissTolerance);
397        Environment::GetSingleton()->GetIntValue("VspTree.Termination.maxViewCells", mMaxViewCells);
398
399        //-- max cost ratio for early tree termination
400        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.maxCostRatio", mTermMaxCostRatio);
401
402        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
403        Environment::GetSingleton()->GetIntValue("VspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
404
405        //-- factors for bsp tree split plane heuristics
406        Environment::GetSingleton()->GetFloatValue("VspTree.Termination.ct_div_ci", mCtDivCi);
407
408        //-- partition criteria
409        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.epsilon", mEpsilon);
410        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
411
412        // if only the driving axis is used for axis aligned split
413        Environment::GetSingleton()->GetBoolValue("VspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
414       
415        Environment::GetSingleton()->GetIntValue("VspTree.maxTests", mMaxTests);
416        Environment::GetSingleton()->GetFloatValue("VspTree.maxStaticMemory", mMaxMemory);
417
418        Environment::GetSingleton()->GetBoolValue("VspTree.useCostHeuristics", mUseCostHeuristics);
419        Environment::GetSingleton()->GetBoolValue("VspTree.simulateOctree", mCirculatingAxis);
420       
421        Environment::GetSingleton()->GetBoolValue("VspTree.useKdPvsForHeuristics", mUseKdPvsForHeuristics);
422
423        char subdivisionStatsLog[100];
424        Environment::GetSingleton()->GetStringValue("VspTree.subdivisionStats", subdivisionStatsLog);
425        mSubdivisionStats.open(subdivisionStatsLog);
426
427        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.minBand", mMinBand);
428        Environment::GetSingleton()->GetFloatValue("VspTree.Construction.maxBand", mMaxBand);
429       
430
431        //-- debug output
432
433        Debug << "******* VSP options ******** " << endl;
434
435    Debug << "max depth: " << mTermMaxDepth << endl;
436        Debug << "min PVS: " << mTermMinPvs << endl;
437        Debug << "min probabiliy: " << mTermMinProbability << endl;
438        Debug << "min rays: " << mTermMinRays << endl;
439        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
440        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
441        Debug << "miss tolerance: " << mTermMissTolerance << endl;
442        Debug << "max view cells: " << mMaxViewCells << endl;
443        Debug << "randomize: " << randomize << endl;
444
445        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
446        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
447        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
448        Debug << "max memory: " << mMaxMemory << endl;
449        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
450        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
451        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
452
453        Debug << "circulating axis: " << mCirculatingAxis << endl;
454        Debug << "minband: " << mMinBand << endl;
455        Debug << "maxband: " << mMaxBand << endl;
456
457        if (!mUseKdPvsForHeuristics)
458                Debug << "pvs heuristics: per object" << endl;
459        else
460                Debug << "pvs heuristics: per kd node" << endl;
461
462        mSplitCandidates = new vector<SortableEntry>;
463
464        Debug << endl;
465}
466
467
468VspViewCell *VspTree::GetOutOfBoundsCell()
469{
470        return mOutOfBoundsCell;
471}
472
473
474VspViewCell *VspTree::GetOrCreateOutOfBoundsCell()
475{
476        if (!mOutOfBoundsCell)
477        {
478                mOutOfBoundsCell = new VspViewCell();
479                mOutOfBoundsCell->SetId(-1);
480                mOutOfBoundsCell->SetValid(false);
481        }
482
483        return mOutOfBoundsCell;
484}
485
486
487const VspTreeStatistics &VspTree::GetStatistics() const
488{
489        return mVspStats;
490}
491
492
493VspTree::~VspTree()
494{
495        DEL_PTR(mRoot);
496        DEL_PTR(mSplitCandidates);
497}
498
499
500void VspTree::ComputeBoundingBox(const VssRayContainer &rays,
501                                                                 AxisAlignedBox3 *forcedBoundingBox)
502{       
503        if (forcedBoundingBox)
504        {
505                mBoundingBox = *forcedBoundingBox;
506        }
507        else // compute vsp tree bounding box
508        {
509                mBoundingBox.Initialize();
510
511                VssRayContainer::const_iterator rit, rit_end = rays.end();
512
513                //-- compute bounding box
514        for (rit = rays.begin(); rit != rit_end; ++ rit)
515                {
516                        VssRay *ray = *rit;
517
518                        // compute bounding box of view space
519                        mBoundingBox.Include(ray->GetTermination());
520                        mBoundingBox.Include(ray->GetOrigin());
521                }
522        }
523}
524
525
526void VspTree::AddSubdivisionStats(const int viewCells,
527                                                                  const float renderCostDecr,
528                                                                  const float totalRenderCost,
529                                                                  const float avgRenderCost)
530{
531        mSubdivisionStats
532                        << "#ViewCells\n" << viewCells << endl
533                        << "#RenderCostDecrease\n" << renderCostDecr << endl
534                        << "#TotalRenderCost\n" << totalRenderCost << endl
535                        << "#AvgRenderCost\n" << avgRenderCost << endl;
536}
537
538
539// TODO: return memory usage in MB
540float VspTree::GetMemUsage() const
541{
542        return (float)
543                 (sizeof(VspTree) +
544                  mVspStats.Leaves() * sizeof(VspLeaf) +
545                  mCreatedViewCells * sizeof(VspViewCell) +
546                  mVspStats.pvs * sizeof(ObjectPvsData) +
547                  mVspStats.Interior() * sizeof(VspInterior) +
548                  mVspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
549}
550
551
552bool VspTree::LocalTerminationCriteriaMet(const VspTraversalData &data) const
553{
554        const bool localTerminationCriteriaMet = (
555                ((int)data.mRays->size() <= mTermMinRays) ||
556                (data.mPvs <= mTermMinPvs)   ||
557                (data.mProbability <= mTermMinProbability) ||
558                (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
559                (data.mDepth >= mTermMaxDepth)
560                );
561
562        if (0 && localTerminationCriteriaMet)
563        {
564                Debug << "********local termination *********" << endl;
565                Debug << "rays: " << (int)data.mRays->size() << "  " << mTermMinRays << endl;
566                Debug << "pvs: " << data.mPvs << " " << mTermMinPvs << endl;
567                Debug << "p: " <<  data.mProbability << " " << mTermMinProbability << endl;
568                Debug << "avg contri: " << data.GetAvgRayContribution() << " " << mTermMaxRayContribution << endl;
569                Debug << "depth " << data.mDepth << " " << mTermMaxDepth << endl;
570        }
571
572        return localTerminationCriteriaMet;
573               
574}
575
576
577bool VspTree::GlobalTerminationCriteriaMet(const VspTraversalData &data) const
578{
579        const bool terminationCriteriaMet = (
580                // mOutOfMemory ||
581                (mVspStats.Leaves() >= mMaxViewCells) ||
582        (mGlobalCostMisses >= mTermGlobalCostMissTolerance)
583                );
584
585        if (0 && terminationCriteriaMet)
586        {
587                Debug << "********* terminationCriteriaMet *********" << endl;
588                Debug << "cost misses: " << mGlobalCostMisses << " " << mTermGlobalCostMissTolerance << endl;
589                Debug << "leaves: " << mVspStats.Leaves() << " " <<  mMaxViewCells << endl;
590        }
591        return terminationCriteriaMet;
592}
593
594
595void VspTree::CreateViewCell(VspTraversalData &tData, const bool updatePvs)
596{
597        //-- create new view cell
598        VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
599       
600        VspViewCell *viewCell = new VspViewCell();
601    leaf->SetViewCell(viewCell);
602       
603        int conSamp = 0;
604        float sampCon = 0.0f;
605
606        //-- update pvs
607        if (updatePvs)
608        {
609                AddSamplesToPvs(leaf, *tData.mRays, sampCon, conSamp);
610
611                // update scalar pvs size value
612                mViewCellsManager->SetScalarPvsSize(viewCell, viewCell->GetPvs().CountPvs());
613
614                mVspStats.contributingSamples += conSamp;
615                mVspStats.sampleContributions += (int)sampCon;
616        }
617
618
619        //-- store additional info
620        if (mStoreRays)
621        {
622                RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
623
624                for (it = tData.mRays->begin(); it != it_end; ++ it)
625                {
626                        (*it).mRay->Ref();                     
627                        leaf->mVssRays.push_back((*it).mRay);
628                }
629        }
630               
631        //
632        viewCell->mLeaf = leaf;
633
634        viewCell->SetVolume(tData.mProbability);
635    leaf->mProbability = tData.mProbability;
636}
637
638
639void VspTree::EvalSubdivisionStats(const VspTraversalData &tData,
640                                                                   const VspTraversalData &tFrontData,
641                                                                   const VspTraversalData &tBackData
642                                                                   )
643{
644        const float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
645        const float cBack = (float)tBackData.mPvs * tBackData.mProbability;
646        const float cData = (float)tData.mPvs * tData.mProbability;
647       
648        const float costDecr =
649                (cFront + cBack - cData) / mBoundingBox.GetVolume();
650
651        mTotalCost += costDecr;
652        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
653
654        AddSubdivisionStats(mVspStats.Leaves(),
655                                                -costDecr,
656                                                mTotalCost,
657                                                (float)mTotalPvsSize / (float)mVspStats.Leaves());
658}
659
660
661VspNode *VspTree::Subdivide(SplitQueue &tQueue,
662                                                        SplitCandidate *splitCandidate,
663                                                        const bool globalCriteriaMet)
664{
665        // doto remove dynamic cast
666        VspSplitCandidate *sc = dynamic_cast<VspSplitCandidate *>(splitCandidate);
667        VspTraversalData &tData = sc->mParentData;
668
669        VspNode *newNode = tData.mNode;
670
671        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
672        {       
673                VspTraversalData tFrontData;
674                VspTraversalData tBackData;
675
676                //-- continue subdivision
677               
678                // create new interior node and two leaf node
679                const AxisAlignedPlane splitPlane = sc->mSplitPlane;
680                newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData);
681       
682                const int maxCostMisses = sc->mMaxCostMisses;
683
684
685                // how often was max cost ratio missed in this branch?
686                tFrontData.mMaxCostMisses = maxCostMisses;
687                tBackData.mMaxCostMisses = maxCostMisses;
688                       
689       
690                if (1) //-- subdivision statistics
691                        EvalSubdivisionStats(tData, tFrontData, tBackData);
692               
693                //-- evaluate new split candidates for global greedy cost heuristics
694
695                VspSplitCandidate *frontCandidate = new VspSplitCandidate(tFrontData);
696                VspSplitCandidate *backCandidate = new VspSplitCandidate(tBackData);
697
698                EvalSplitCandidate(*frontCandidate);
699                EvalSplitCandidate(*backCandidate);
700
701                tQueue.Push(frontCandidate);
702                tQueue.Push(backCandidate);
703               
704                // delete old view cell
705                delete tData.mNode->mViewCell;
706
707                // delete old leaf node
708                DEL_PTR(tData.mNode);
709        }
710
711        if (newNode->IsLeaf()) // subdivision terminated
712        {
713                //CreateViewCell(tData); // create new view cell
714
715                VspLeaf *leaf = dynamic_cast<VspLeaf *>(newNode);
716                ViewCell *viewCell = leaf->GetViewCell();
717
718                int conSamp = 0;
719                float sampCon = 0.0f;
720
721
722                //-- store pvs optained from rays
723                AddSamplesToPvs(leaf, *tData.mRays, sampCon, conSamp);
724
725                // update scalar pvs size value
726                mViewCellsManager->SetScalarPvsSize(viewCell, viewCell->GetPvs().CountPvs());
727
728                mVspStats.contributingSamples += conSamp;
729                mVspStats.sampleContributions += (int)sampCon;
730
731                //-- store additional info
732                if (mStoreRays)
733                {
734                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
735                        for (it = tData.mRays->begin(); it != it_end; ++ it)
736                        {
737                                (*it).mRay->Ref();                     
738                                dynamic_cast<VspLeaf *>(newNode)->mVssRays.push_back((*it).mRay);
739                        }
740                }
741
742                // finally evaluate statistics for this leaf
743                EvaluateLeafStats(tData);
744        }
745
746        //-- cleanup
747        tData.Clear();
748       
749        return newNode;
750}
751
752
753void VspTree::EvalSplitCandidate(VspSplitCandidate &splitCandidate)
754{
755        float frontProb;
756        float backProb;
757       
758        VspLeaf *leaf = dynamic_cast<VspLeaf *>(splitCandidate.mParentData.mNode);
759       
760        // compute locally best split plane
761        const bool success =
762                SelectSplitPlane(splitCandidate.mParentData, splitCandidate.mSplitPlane, frontProb, backProb);
763
764        float oldRenderCost;
765
766        // compute global decrease in render cost
767        const float renderCostDecr = EvalRenderCostDecrease(splitCandidate.mSplitPlane,
768                                                                                                                splitCandidate.mParentData,
769                                                                                                                oldRenderCost);
770
771        splitCandidate.SetRenderCostDecrease(renderCostDecr);
772
773#if 0
774        const float priority = (float)-data.mDepth;
775#else   
776
777        // take render cost of node into account
778        // otherwise danger of being stuck in a local minimum!!
779        const float factor = mRenderCostDecreaseWeight;
780        const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;
781#endif
782       
783        splitCandidate.SetPriority(priority);
784
785        // max cost threshold violated?
786        splitCandidate.mMaxCostMisses =
787                success ? splitCandidate.mParentData.mMaxCostMisses : splitCandidate.mParentData.mMaxCostMisses + 1;
788        //Debug << "p: " << tData.mNode << " depth: " << tData.mDepth << endl;
789}
790
791
792VspInterior *VspTree::SubdivideNode(const AxisAlignedPlane &splitPlane,
793                                                                        VspTraversalData &tData,
794                                                                        VspTraversalData &frontData,
795                                                                        VspTraversalData &backData)
796{
797        VspLeaf *leaf = dynamic_cast<VspLeaf *>(tData.mNode);
798       
799        //-- the front and back traversal data is filled with the new values
800
801        frontData.mDepth = tData.mDepth + 1;
802        backData.mDepth = tData.mDepth + 1;
803
804        frontData.mRays = new RayInfoContainer();
805        backData.mRays = new RayInfoContainer();
806
807        //-- subdivide rays
808        SplitRays(splitPlane,
809                          *tData.mRays,
810                          *frontData.mRays,
811                          *backData.mRays);
812       
813        //Debug << "f: " << frontData.mRays->size() << " b: " << backData.mRays->size() << "d: " << tData.mRays->size() << endl;
814        //-- compute pvs
815        frontData.mPvs = ComputePvsSize(*frontData.mRays);
816        backData.mPvs = ComputePvsSize(*backData.mRays);
817
818        // split front and back node geometry and compute area
819        tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition,
820                                                         frontData.mBoundingBox, backData.mBoundingBox);
821
822
823        frontData.mProbability = frontData.mBoundingBox.GetVolume();
824        backData.mProbability = tData.mProbability - frontData.mProbability;
825
826       
827    ///////////////////////////////////////////
828        // subdivide further
829
830        // store maximal and minimal depth
831        if (tData.mDepth > mVspStats.maxDepth)
832        {
833                Debug << "max depth increases to " << tData.mDepth << " at " << mVspStats.Leaves() << " leaves" << endl;
834                mVspStats.maxDepth = tData.mDepth;
835        }
836
837        // two more leaves
838        mVspStats.nodes += 2;
839   
840        VspInterior *interior = new VspInterior(splitPlane);
841
842#ifdef _DEBUG
843        Debug << interior << endl;
844#endif
845
846
847        //-- create front and back leaf
848
849        VspInterior *parent = leaf->GetParent();
850
851        // replace a link from node's parent
852        if (parent)
853        {
854                parent->ReplaceChildLink(leaf, interior);
855                interior->SetParent(parent);
856
857                // remove "parent" view cell from pvs of all objects (traverse trough rays)
858                RemoveParentViewCellReferences(tData.mNode->GetViewCell());
859        }
860        else // new root
861        {
862                mRoot = interior;
863        }
864
865        VspLeaf *frontLeaf = new VspLeaf(interior);
866        VspLeaf *backLeaf = new VspLeaf(interior);
867
868        // and setup child links
869        interior->SetupChildLinks(frontLeaf, backLeaf);
870       
871        // add bounding box
872        interior->SetBoundingBox(tData.mBoundingBox);
873
874        // set front and back leaf
875        frontData.mNode = frontLeaf;
876        backData.mNode = backLeaf;
877
878        // explicitely create front and back view cell
879        CreateViewCell(frontData, false);
880        CreateViewCell(backData, false);
881
882       
883        // create front and back view cell
884        // add front and back view cell to "Potentially Visbilie View Cells"
885        // of the objects in front and back pvs
886
887        //AddViewCellReferences(frontLeaf->GetViewCell());
888        //AddViewCellReferences(backLeaf->GetViewCell());
889
890        interior->mTimeStamp = mTimeStamp ++;
891       
892        return interior;
893}
894
895
896KdNode *VspTree::GetKdNode(VssRay &ray, const bool termination) const
897{
898        if (termination)
899        {
900                if (!ray.mTerminationNode)
901                {
902                        ray.mTerminationNode = mOspTree->GetRoot();
903                }
904
905                ray.mTerminationNode =
906                        mOspTree->GetLeaf(ray.mTermination, ray.mTerminationNode);
907
908                return ray.mTerminationNode;
909        }
910        else
911        {
912                if (!ray.mOriginNode)
913                {
914                        ray.mOriginNode = mOspTree->GetRoot();
915                }
916
917                ray.mOriginNode =
918                        mOspTree->GetLeaf(ray.mOrigin, ray.mOriginNode);
919
920                return ray.mOriginNode;
921        }
922}
923
924
925void VspTree::RemoveParentViewCellReferences(ViewCell *parent) const
926{
927        KdLeaf::NewMail();
928
929        // remove the parents from the object pvss
930        ObjectPvsMap::const_iterator oit, oit_end = parent->GetPvs().mEntries.end();
931
932        for (oit = parent->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
933        {
934                Intersectable *object = (*oit).first;
935                // HACK: make sure that the view cell is removed from the pvs
936                const float high_contri = 9999999;
937
938                // remove reference count of view cells
939                object->mViewCellPvs.RemoveSample(parent, high_contri);
940        }
941}
942
943
944void VspTree::AddViewCellReferences(ViewCell *vc) const
945{
946        KdLeaf::NewMail();
947
948        // Add front view cell to the object pvsss
949        ObjectPvsMap::const_iterator oit, oit_end = vc->GetPvs().mEntries.end();
950
951        for (oit = vc->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
952        {
953                Intersectable *object = (*oit).first;
954
955                // increase reference count of view cells
956                object->mViewCellPvs.AddSample(vc, 1);
957        }
958}
959
960
961bool VspTree::AddKdLeafToPvs(KdLeaf *leaf,
962                                                         ViewCell *vc,
963                                                         float &pdf,
964                                                         float &contribution)
965{
966        bool contri = false;
967
968#if 1 // add kd intersecable to pvs
969        KdIntersectable *kdObj = mOspTree->GetOrCreateKdIntersectable(leaf);
970       
971        if (vc->AddPvsSample(kdObj, pdf, contribution))
972        {
973                return true;
974        }
975
976#else // add all objects of kd node
977
978        pdf = 0;
979        contribution = 0;
980
981        ObjectContainer::const_iterator it, it_end = leaf->mObjects.end();
982
983        for (it = leaf->mObjects.begin(); it != it_end; ++ it)
984        {
985                Intersectable *object = *it;                                           
986
987                float newpdf;
988                float newcontri;
989                                               
990                if (vc->AddPvsSample(object, newpdf, newcontri))
991                {
992                        contri = true;
993                }
994
995                pdf += newPdf;
996                newContri += contribution;
997        }
998
999#endif
1000
1001        return contri;
1002}
1003
1004void VspTree::AddSamplesToPvs(VspLeaf *leaf,
1005                                                         const RayInfoContainer &rays,
1006                                                         float &sampleContributions,
1007                                                         int &contributingSamples)
1008{
1009        sampleContributions = 0;
1010        contributingSamples = 0;
1011 
1012        RayInfoContainer::const_iterator it, it_end = rays.end();
1013 
1014        ViewCellLeaf *vc = leaf->GetViewCell();
1015 
1016        // add contributions from samples to the PVS
1017        for (it = rays.begin(); it != it_end; ++ it)
1018        {
1019                float sc = 0.0f;
1020                VssRay *ray = (*it).mRay;
1021
1022                bool madeContrib = false;
1023                float contribution;
1024
1025                Intersectable *obj = ray->mTerminationObject;
1026
1027                if (obj)
1028                {
1029                        if (mStoreKdPvs)
1030                        {
1031                                // potentially visible kd cells
1032                                KdLeaf *leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
1033                                AddKdLeafToPvs(leaf, vc, ray->mPdf, contribution);
1034                        }
1035                        else
1036                        {
1037                                if (vc->AddPvsSample(obj, ray->mPdf, contribution))
1038                                {
1039                                        madeContrib = true;
1040                                }
1041                        }
1042
1043                        sc += contribution;
1044                }
1045
1046                obj = ray->mOriginObject;
1047
1048                if (obj)
1049                {
1050                        // potentially visible kd cells
1051                        if (mUseKdPvsForHeuristics)
1052                        {
1053                                KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);       
1054                                AddKdLeafToPvs(leaf, vc, ray->mPdf, contribution);
1055                        }
1056                        else
1057                        {
1058                                if (vc->AddPvsSample(obj, ray->mPdf, contribution))
1059                                {
1060                                        madeContrib = true;
1061                                }
1062                        }
1063
1064                        sc += contribution;
1065                }
1066
1067                if (madeContrib)
1068                {
1069                        ++ contributingSamples;
1070                }
1071
1072                // store rays for visualization
1073                if (0) leaf->mVssRays.push_back(new VssRay(*ray));
1074        }
1075}
1076
1077
1078void VspTree::SortSplitCandidates(const RayInfoContainer &rays,
1079                                                                  const int axis,
1080                                                                  float minBand,
1081                                                                  float maxBand)
1082{
1083        mSplitCandidates->clear();
1084
1085        int requestedSize = 2 * (int)(rays.size());
1086
1087        // creates a sorted split candidates array
1088        if (mSplitCandidates->capacity() > 500000 &&
1089                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
1090        {
1091        delete mSplitCandidates;
1092                mSplitCandidates = new vector<SortableEntry>;
1093        }
1094
1095        mSplitCandidates->reserve(requestedSize);
1096
1097        float pos;
1098
1099        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1100
1101        //-- insert all queries
1102        for (rit = rays.begin(); rit != rit_end; ++ rit)
1103        {
1104                const bool positive = (*rit).mRay->HasPosDir(axis);
1105                               
1106                pos = (*rit).ExtrapOrigin(axis);
1107               
1108                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
1109                                                                        pos, (*rit).mRay));
1110
1111                pos = (*rit).ExtrapTermination(axis);
1112
1113                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
1114                                                                        pos, (*rit).mRay));
1115        }
1116
1117        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
1118}
1119
1120
1121int VspTree::GetPvsContribution(Intersectable *object) const
1122{
1123    int pvsContri = 0;
1124
1125        KdPvsMap::const_iterator kit, kit_end = object->mKdPvs.mEntries.end();
1126
1127        Intersectable::NewMail();
1128
1129        // Search kd leaves this object is attached to
1130        for (kit = object->mKdPvs.mEntries.begin(); kit != kit_end; ++ kit)
1131        {
1132                KdNode *l = (*kit).first;
1133
1134                // new object found during sweep
1135                // => increase pvs contribution of this kd node
1136                if (!l->Mailed())
1137                {
1138                        l->Mail();
1139                        ++ pvsContri;
1140                }
1141        }
1142
1143        return pvsContri;
1144}
1145
1146
1147int VspTree::PrepareHeuristics(KdLeaf *leaf)
1148{       
1149        int pvsSize = 0;
1150       
1151        if (!leaf->Mailed())
1152        {
1153                leaf->Mail();
1154                leaf->mCounter = 1;
1155
1156                // add objects without the objects which are in several kd leaves
1157                pvsSize += (int)(leaf->mObjects.size() - leaf->mMultipleObjects.size());
1158                //Debug << "here11 " << (int)leaf->mObjects.size() << " " << leaf->mMultipleObjects.size() << endl;
1159        }
1160        else
1161        {
1162                ++ leaf->mCounter;
1163        }
1164
1165        //-- the objects belonging to several leaves must be handled seperately
1166
1167        ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1168
1169        for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1170        {
1171                Intersectable *object = *oit;
1172                                               
1173                if (!object->Mailed())
1174                {
1175                        object->Mail();
1176                        object->mCounter = 1;
1177
1178                        ++ pvsSize;
1179                }
1180                else
1181                {
1182                        ++ object->mCounter;
1183                }
1184        }
1185       
1186        return pvsSize;
1187}
1188
1189
1190int VspTree::PrepareHeuristics(const RayInfoContainer &rays)
1191{       
1192        Intersectable::NewMail();
1193        KdNode::NewMail();
1194
1195        int pvsSize = 0;
1196
1197        RayInfoContainer::const_iterator ri, ri_end = rays.end();
1198
1199        //-- set all kd nodes / objects as belonging to the front pvs
1200
1201        for (ri = rays.begin(); ri != ri_end; ++ ri)
1202        {
1203                VssRay *ray = (*ri).mRay;
1204                Intersectable *oObject = ray->mOriginObject;
1205
1206                if (oObject)
1207                {
1208                        if (!mUseKdPvsForHeuristics)
1209                        {
1210                                if (!oObject->Mailed())
1211                                {
1212                                        oObject->Mail();
1213                                        oObject->mCounter = 1;
1214                                        ++ pvsSize;
1215                                }
1216                                else
1217                                {
1218                                        ++ oObject->mCounter;
1219                                }
1220                        }
1221                        else
1222                        {
1223                                KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);
1224                                pvsSize += PrepareHeuristics(leaf);     
1225                        }       
1226                }
1227
1228                Intersectable *tObject = (*ri).mRay->mTerminationObject;
1229
1230                if (tObject)
1231                {
1232                        if (!mUseKdPvsForHeuristics)
1233                        {
1234                                if (!tObject->Mailed())
1235                                {
1236                                        tObject->Mail();
1237                                        tObject->mCounter = 1;
1238                                        ++ pvsSize;
1239                                }
1240                                else
1241                                {
1242                                        ++ tObject->mCounter;
1243                                }
1244                        }
1245                        else
1246                        {
1247                                KdLeaf *leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
1248                                pvsSize += PrepareHeuristics(leaf);     
1249                        }       
1250                }
1251        }
1252
1253        return pvsSize;
1254}
1255
1256
1257void VspTree::RemoveContriFromPvs(KdLeaf *leaf, int &pvs) const
1258{
1259        // leaf falls out of right pvs
1260        if (-- leaf->mCounter == 0)
1261        {
1262                pvs -= ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
1263        }
1264
1265        //-- handle objects which are in several kd leaves separately
1266        ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1267
1268        for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1269        {
1270                Intersectable *object = *oit;
1271
1272                if (-- object->mCounter == 0)
1273                {
1274                        -- pvs;
1275                }
1276        }
1277}
1278
1279
1280void VspTree::AddContriToPvs(KdLeaf *leaf, int &pvs) const
1281{
1282        if (leaf->Mailed())
1283                return;
1284       
1285        leaf->Mail();
1286
1287        // add objects without those which are part of several kd leaves
1288        pvs += ((int)leaf->mObjects.size() - (int)leaf->mMultipleObjects.size());
1289
1290        //-- handle objects of several kd leaves separately
1291        ObjectContainer::const_iterator oit, oit_end = leaf->mMultipleObjects.end();
1292
1293        for (oit = leaf->mMultipleObjects.begin(); oit != oit_end; ++ oit)
1294        {
1295                Intersectable *object = *oit;
1296
1297                // object not previously in pvs
1298                if (!object->Mailed())
1299                {
1300                        object->Mail();
1301                        ++ pvs;
1302                }
1303        }                                       
1304}
1305
1306
1307void VspTree::EvalPvsIncr(const SortableEntry &ci,
1308                                                  int &pvsLeft,
1309                                                  int &pvsRight) const
1310{
1311        VssRay *ray = ci.ray;
1312
1313        Intersectable *oObject = ray->mOriginObject;
1314        Intersectable *tObject = ray->mTerminationObject;
1315
1316        if (oObject)
1317        {       
1318                if (!mUseKdPvsForHeuristics)
1319                {
1320                        if (ci.type == SortableEntry::ERayMin)
1321                        {
1322                                if (!oObject->Mailed())
1323                                {
1324                                        oObject->Mail();
1325                                        ++ pvsLeft;
1326                                }
1327                        }
1328                        else if (ci.type == SortableEntry::ERayMax)
1329                        {
1330                                if (-- oObject->mCounter == 0)
1331                                        -- pvsRight;
1332                        }
1333                }
1334                else // per kd node
1335                {
1336                        KdLeaf *node = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);
1337
1338                        // add contributions of the kd nodes
1339                        if (ci.type == SortableEntry::ERayMin)
1340                        {
1341                                AddContriToPvs(node, pvsLeft);
1342                        }
1343                        else if (ci.type == SortableEntry::ERayMax)
1344                        {
1345                                RemoveContriFromPvs(node, pvsRight);
1346                        }
1347                }
1348                       
1349
1350        }
1351       
1352        if (tObject)
1353        {       
1354                if (!mUseKdPvsForHeuristics)
1355                {
1356                        if (ci.type == SortableEntry::ERayMin)
1357                        {
1358                                if (!tObject->Mailed())
1359                                {
1360                                        tObject->Mail();
1361                                        ++ pvsLeft;
1362                                }
1363                        }
1364                        else if (ci.type == SortableEntry::ERayMax)
1365                        {
1366                                if (-- tObject->mCounter == 0)
1367                                        -- pvsRight;
1368                        }
1369                }
1370                else // per kd node
1371                {
1372                        KdLeaf *node = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
1373
1374                        if (ci.type == SortableEntry::ERayMin)
1375                        {
1376                                AddContriToPvs(node, pvsLeft);
1377                        }
1378                        else if (ci.type == SortableEntry::ERayMax)
1379                        {
1380                                RemoveContriFromPvs(node, pvsRight);
1381                        }
1382                }
1383        }
1384}
1385
1386
1387float VspTree::EvalLocalCostHeuristics(const VspTraversalData &tData,
1388                                                                           const AxisAlignedBox3 &box,
1389                                                                           const int axis,
1390                                                                           float &position)
1391{
1392        RayInfoContainer usedRays;
1393
1394        if (mMaxTests < (int)tData.mRays->size())
1395        {
1396                GetRayInfoSets(*tData.mRays, mMaxTests, usedRays);
1397        }
1398        else
1399        {
1400                usedRays = *tData.mRays;
1401        }
1402
1403        int pvsSize = tData.mPvs;
1404
1405        const float minBox = box.Min(axis);
1406        const float maxBox = box.Max(axis);
1407
1408        const float sizeBox = maxBox - minBox;
1409
1410        const float minBand = minBox + mMinBand * sizeBox;
1411        const float maxBand = minBox + mMaxBand * sizeBox;
1412
1413        SortSplitCandidates(usedRays, axis, minBand, maxBand);
1414
1415        // prepare the sweep
1416        // (note: returns pvs size, so there would be no need
1417        // to give pvs size as argument)
1418        pvsSize = PrepareHeuristics(usedRays);
1419
1420        // go through the lists, count the number of objects left and right
1421        // and evaluate the following cost funcion:
1422        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
1423
1424        int pvsl = 0;
1425        int pvsr = pvsSize;
1426
1427        int pvsBack = pvsl;
1428        int pvsFront = pvsr;
1429
1430        float sum = (float)pvsSize * sizeBox;
1431        float minSum = 1e20f;
1432
1433       
1434        // if no good split can be found, take mid split
1435        position = minBox + 0.5f * sizeBox;
1436       
1437        // the relative cost ratio
1438        float ratio = 99999999.0f;
1439        bool splitPlaneFound = false;
1440
1441        Intersectable::NewMail();
1442        KdLeaf::NewMail();
1443
1444
1445        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
1446
1447        //-- traverse through visibility events
1448
1449        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
1450        {
1451                EvalPvsIncr(*ci, pvsl, pvsr);
1452
1453                // Note: sufficient to compare size of bounding boxes of front and back side?
1454                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
1455                {
1456
1457                        float currentPos;
1458                       
1459                        // HACK: current positition is BETWEEN visibility events
1460                        if (0 && ((ci + 1) != ci_end))
1461                                currentPos = ((*ci).value + (*(ci + 1)).value) * 0.5f;
1462                        else
1463                                currentPos = (*ci).value;                       
1464
1465                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
1466                        //Debug  << "pos=" << (*ci).value << "\t newpos=" << currentPos << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << "\t cost= " << sum << endl;
1467
1468                        if (sum < minSum)
1469                        {
1470                                splitPlaneFound = true;
1471
1472                                minSum = sum;
1473                                position = currentPos;
1474                               
1475                                pvsBack = pvsl;
1476                                pvsFront = pvsr;
1477                        }
1478                }
1479        }
1480       
1481       
1482        //-- compute cost
1483
1484        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1485        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1486
1487        const float pOverall = sizeBox;
1488        const float pBack = position - minBox;
1489        const float pFront = maxBox - position;
1490
1491        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
1492    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1493        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1494       
1495        const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
1496        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1497
1498        if (splitPlaneFound)
1499        {
1500                ratio = newRenderCost / oldRenderCost;
1501        }
1502       
1503        //if (axis != 1)
1504        //Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1505        //      <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
1506
1507        const float volRatio = tData.mBoundingBox.GetVolume() / (sizeBox * mBoundingBox.GetVolume());
1508
1509        Debug << "\n§§§§ eval local cost §§§§" << endl
1510                  << "back pvs: " << penaltyBack << " front pvs: " << penaltyFront << " total pvs: " << penaltyOld << endl
1511                  << "back p: " << pBack * volRatio << " front p " << pFront * volRatio << " p: " << pOverall * volRatio << endl
1512                  << "old rc: " << oldRenderCost * volRatio << " new rc: " << newRenderCost * volRatio << endl
1513                  << "render cost decrease: " << oldRenderCost * volRatio - newRenderCost * volRatio << endl;
1514
1515        return ratio;
1516}
1517
1518
1519float VspTree::SelectSplitPlane(const VspTraversalData &tData,
1520                                                                AxisAlignedPlane &plane,
1521                                                                float &pFront,
1522                                                                float &pBack)
1523{
1524        float nPosition[3];
1525        float nCostRatio[3];
1526        float nProbFront[3];
1527        float nProbBack[3];
1528
1529        // create bounding box of node geometry
1530        AxisAlignedBox3 box = tData.mBoundingBox;
1531               
1532        int sAxis = 0;
1533        int bestAxis = -1;
1534
1535        // if we use some kind of specialised fixed axis
1536    const bool useSpecialAxis =
1537                mOnlyDrivingAxis || mCirculatingAxis;
1538        //Debug << "data: " << tData.mBoundingBox << " pvs " << tData.mPvs << endl;
1539        if (mCirculatingAxis)
1540        {
1541                int parentAxis = 0;
1542                VspNode *parent = tData.mNode->GetParent();
1543
1544                if (parent)
1545                        parentAxis = dynamic_cast<VspInterior *>(parent)->GetAxis();
1546
1547                sAxis = (parentAxis + 1) % 3;
1548        }
1549        else if (mOnlyDrivingAxis)
1550        {
1551                sAxis = box.Size().DrivingAxis();
1552        }
1553       
1554        for (int axis = 0; axis < 3; ++ axis)
1555        {
1556                if (!useSpecialAxis || (axis == sAxis))
1557                {
1558                        if (mUseCostHeuristics)
1559                        {
1560                                //-- place split plane using heuristics
1561                                nCostRatio[axis] =
1562                                        EvalLocalCostHeuristics(tData,
1563                                                                                        box,
1564                                                                                        axis,
1565                                                                                        nPosition[axis]);                       
1566                        }
1567                        else
1568                        {
1569                                //-- split plane position is spatial median
1570                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1571                                nCostRatio[axis] = EvalLocalSplitCost(tData,
1572                                                                                                          box,
1573                                                                                                          axis,
1574                                                                                                          nPosition[axis],
1575                                                                                                          nProbFront[axis],
1576                                                                                                          nProbBack[axis]);
1577                        }
1578                                               
1579                        if (bestAxis == -1)
1580                        {
1581                                bestAxis = axis;
1582                        }
1583                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1584                        {
1585                                bestAxis = axis;
1586                        }
1587                }
1588        }
1589
1590
1591        //-- assign values of best split
1592       
1593        plane.mAxis = bestAxis;
1594        plane.mPosition = nPosition[bestAxis]; // split plane position
1595
1596        pFront = nProbFront[bestAxis];
1597        pBack = nProbBack[bestAxis];
1598
1599        //Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
1600        return nCostRatio[bestAxis];
1601}
1602
1603
1604float VspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
1605                                                                          const VspTraversalData &data,
1606                                                                          float &normalizedOldRenderCost) const
1607{
1608        float pvsFront = 0;
1609        float pvsBack = 0;
1610        float totalPvs = 0;
1611
1612        // probability that view point lies in back / front node
1613        float pOverall = data.mProbability;
1614        float pFront = 0;
1615        float pBack = 0;
1616
1617        const float viewSpaceVol = mBoundingBox.GetVolume();
1618
1619        // create unique ids for pvs heuristics
1620        Intersectable::NewMail(3);
1621       
1622        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1623
1624        for (rit = data.mRays->begin(); rit != rit_end; ++ rit)
1625        {
1626                RayInfo rayInf = *rit;
1627
1628                float t;
1629                VssRay *ray = rayInf.mRay;
1630
1631                const int cf =
1632                        rayInf.ComputeRayIntersection(candidatePlane.mAxis,
1633                                                                                  candidatePlane.mPosition, t);
1634
1635                // find front and back pvs for origing and termination object
1636                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1637                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1638        }
1639
1640
1641        AxisAlignedBox3 frontBox;
1642        AxisAlignedBox3 backBox;
1643
1644        data.mBoundingBox.Split(candidatePlane.mAxis, candidatePlane.mPosition, frontBox, backBox);
1645
1646        pFront = frontBox.GetVolume();
1647        pBack = pOverall - pFront;
1648               
1649
1650        //-- pvs rendering heuristics
1651        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1652        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1653
1654        //-- only render cost heuristics or combined with standard deviation
1655        const float penaltyOld = EvalPvsPenalty((int)totalPvs, lowerPvsLimit, upperPvsLimit);
1656    const float penaltyFront = EvalPvsPenalty((int)pvsFront, lowerPvsLimit, upperPvsLimit);
1657        const float penaltyBack = EvalPvsPenalty((int)pvsBack, lowerPvsLimit, upperPvsLimit);
1658                       
1659        const float oldRenderCost = pOverall * penaltyOld;
1660        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1661
1662        normalizedOldRenderCost = oldRenderCost / viewSpaceVol;
1663
1664        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
1665        const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;
1666       
1667       
1668        Debug << "\n==== eval render cost decrease ===" << endl
1669                  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl
1670                  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol << " p: " << pOverall / viewSpaceVol << endl
1671                  << "old rc: " << normalizedOldRenderCost << " new rc: " << newRenderCost / viewSpaceVol << endl
1672                  << "render cost decrease: " << renderCostDecrease << endl;
1673
1674        return renderCostDecrease;
1675}
1676
1677
1678float VspTree::EvalLocalSplitCost(const VspTraversalData &data,
1679                                                                  const AxisAlignedBox3 &box,
1680                                                                  const int axis,
1681                                                                  const float &position,
1682                                                                  float &pFront,
1683                                                                  float &pBack) const
1684{
1685        float pvsTotal = 0;
1686        float pvsFront = 0;
1687        float pvsBack = 0;
1688       
1689        // create unique ids for pvs heuristics
1690        Intersectable::NewMail(3);
1691
1692        const int pvsSize = data.mPvs;
1693
1694        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1695
1696        // this is the main ray classification loop!
1697        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
1698        {
1699                // determine the side of this ray with respect to the plane
1700                float t;
1701                const int side = (*rit).ComputeRayIntersection(axis, position, t);
1702       
1703                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
1704                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
1705        }
1706
1707        //-- cost heuristics
1708        float pOverall;
1709       
1710        pOverall = data.mProbability;
1711
1712        // we use spatial mid split => simplified computation
1713        pBack = pFront = pOverall * 0.5f;
1714       
1715        const float newCost = pvsBack * pBack + pvsFront * pFront;
1716        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
1717       
1718#ifdef _DEBUG
1719        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
1720        Debug << pFront << " " << pBack << " " << pOverall << endl;
1721#endif
1722
1723        return  (mCtDivCi + newCost) / oldCost;
1724}
1725
1726
1727void VspTree::AddObjToPvs(Intersectable *obj,
1728                                                  const int cf,
1729                                                  float &frontPvs,
1730                                                  float &backPvs,
1731                                                  float &totalPvs) const
1732{
1733        if (!obj)
1734                return;
1735
1736        //const float renderCost = mViewCellsManager->SimpleRay &raynderCost(obj);
1737        const int renderCost = 1;
1738
1739        // object in no pvs => new
1740        if (!obj->Mailed() && !obj->Mailed(1) && !obj->Mailed(2))
1741        {
1742                totalPvs += renderCost;
1743        }
1744
1745        // QUESTION matt: is it safe to assume that
1746        // the object belongs to no pvs in this case?
1747        //if (cf == Ray::COINCIDENT) return;
1748
1749        if (cf >= 0) // front pvs
1750        {
1751                if (!obj->Mailed() && !obj->Mailed(2))
1752                {
1753                        frontPvs += renderCost;
1754               
1755                        // already in back pvs => in both pvss
1756                        if (obj->Mailed(1))
1757                                obj->Mail(2);
1758                        else
1759                                obj->Mail();
1760                }
1761        }
1762
1763        if (cf <= 0) // back pvs
1764        {
1765                if (!obj->Mailed(1) && !obj->Mailed(2))
1766                {
1767                        backPvs += renderCost;
1768               
1769                        // already in front pvs => in both pvss
1770                        if (obj->Mailed())
1771                                obj->Mail(2);
1772                        else
1773                                obj->Mail(1);
1774                }
1775        }
1776}
1777
1778
1779void VspTree::CollectLeaves(vector<VspLeaf *> &leaves,
1780                                                           const bool onlyUnmailed,
1781                                                           const int maxPvsSize) const
1782{
1783        stack<VspNode *> nodeStack;
1784        nodeStack.push(mRoot);
1785
1786        while (!nodeStack.empty())
1787        {
1788                VspNode *node = nodeStack.top();
1789                nodeStack.pop();
1790               
1791                if (node->IsLeaf())
1792                {
1793                        // test if this leaf is in valid view space
1794                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1795                        if (leaf->TreeValid() &&
1796                                (!onlyUnmailed || !leaf->Mailed()) &&
1797                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().CountPvs() <= maxPvsSize)))
1798                        {
1799                                leaves.push_back(leaf);
1800                        }
1801                }
1802                else
1803                {
1804                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1805
1806                        nodeStack.push(interior->GetBack());
1807                        nodeStack.push(interior->GetFront());
1808                }
1809        }
1810}
1811
1812
1813AxisAlignedBox3 VspTree::GetBoundingBox() const
1814{
1815        return mBoundingBox;
1816}
1817
1818
1819VspNode *VspTree::GetRoot() const
1820{
1821        return mRoot;
1822}
1823
1824
1825void VspTree::EvaluateLeafStats(const VspTraversalData &data)
1826{
1827        // the node became a leaf -> evaluate stats for leafs
1828        VspLeaf *leaf = dynamic_cast<VspLeaf *>(data.mNode);
1829
1830
1831        if (data.mPvs > mVspStats.maxPvs)
1832        {
1833                mVspStats.maxPvs = data.mPvs;
1834        }
1835
1836        mVspStats.pvs += data.mPvs;
1837
1838        if (data.mDepth < mVspStats.minDepth)
1839        {
1840                mVspStats.minDepth = data.mDepth;
1841        }
1842       
1843        if (data.mDepth >= mTermMaxDepth)
1844        {
1845        ++ mVspStats.maxDepthNodes;
1846                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
1847        }
1848
1849        // accumulate rays to compute rays /  leaf
1850        mVspStats.accumRays += (int)data.mRays->size();
1851
1852        if (data.mPvs < mTermMinPvs)
1853                ++ mVspStats.minPvsNodes;
1854
1855        if ((int)data.mRays->size() < mTermMinRays)
1856                ++ mVspStats.minRaysNodes;
1857
1858        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1859                ++ mVspStats.maxRayContribNodes;
1860
1861        if (data.mProbability <= mTermMinProbability)
1862                ++ mVspStats.minProbabilityNodes;
1863       
1864        // accumulate depth to compute average depth
1865        mVspStats.accumDepth += data.mDepth;
1866
1867        ++ mCreatedViewCells;
1868
1869#ifdef _DEBUG
1870        Debug << "BSP stats: "
1871                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1872                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1873                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1874                  << "#pvs: " << leaf->GetViewCell()->GetPvs().CountPvs() << "), "
1875                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1876#endif
1877}
1878
1879
1880void VspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
1881{
1882        ViewCell::NewMail();
1883        CollectViewCells(mRoot, onlyValid, viewCells, true);
1884}
1885
1886
1887void VspTree::CollapseViewCells()
1888{
1889// TODO matt
1890#if HAS_TO_BE_REDONE
1891        stack<VspNode *> nodeStack;
1892
1893        if (!mRoot)
1894                return;
1895
1896        nodeStack.push(mRoot);
1897       
1898        while (!nodeStack.empty())
1899        {
1900                VspNode *node = nodeStack.top();
1901                nodeStack.pop();
1902               
1903                if (node->IsLeaf())
1904        {
1905                        BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1906
1907                        if (!viewCell->GetValid())
1908                        {
1909                                BspViewCell *viewCell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
1910       
1911                                ViewCellContainer leaves;
1912                                mViewCellsTree->CollectLeaves(viewCell, leaves);
1913
1914                                ViewCellContainer::const_iterator it, it_end = leaves.end();
1915
1916                                for (it = leaves.begin(); it != it_end; ++ it)
1917                                {
1918                                        VspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
1919                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
1920                                        ++ mVspStats.invalidLeaves;
1921                                }
1922
1923                                // add to unbounded view cell
1924                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
1925                                DEL_PTR(viewCell);
1926                        }
1927                }
1928                else
1929                {
1930                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1931               
1932                        nodeStack.push(interior->GetFront());
1933                        nodeStack.push(interior->GetBack());
1934                }
1935        }
1936
1937        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
1938#endif
1939}
1940
1941
1942void VspTree::CollectRays(VssRayContainer &rays)
1943{
1944        vector<VspLeaf *> leaves;
1945
1946        vector<VspLeaf *>::const_iterator lit, lit_end = leaves.end();
1947
1948        for (lit = leaves.begin(); lit != lit_end; ++ lit)
1949        {
1950                VspLeaf *leaf = *lit;
1951                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
1952
1953                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
1954                        rays.push_back(*rit);
1955        }
1956}
1957
1958
1959void VspTree::SetViewCellsManager(ViewCellsManager *vcm)
1960{
1961        mViewCellsManager = vcm;
1962}
1963
1964
1965void VspTree::ValidateTree()
1966{
1967        mVspStats.invalidLeaves = 0;
1968
1969        stack<VspNode *> nodeStack;
1970
1971        if (!mRoot)
1972                return;
1973
1974        nodeStack.push(mRoot);
1975
1976        while (!nodeStack.empty())
1977        {
1978                VspNode *node = nodeStack.top();
1979                nodeStack.pop();
1980               
1981                if (node->IsLeaf())
1982                {
1983                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
1984
1985                        if (!leaf->GetViewCell()->GetValid())
1986                                ++ mVspStats.invalidLeaves;
1987
1988                        // validity flags don't match => repair
1989                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
1990                        {
1991                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
1992                                PropagateUpValidity(leaf);
1993                        }
1994                }
1995                else
1996                {
1997                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
1998               
1999                        nodeStack.push(interior->GetFront());
2000                        nodeStack.push(interior->GetBack());
2001                }
2002        }
2003
2004        Debug << "invalid leaves: " << mVspStats.invalidLeaves << endl;
2005}
2006
2007
2008
2009void VspTree::CollectViewCells(VspNode *root,
2010                                                                  bool onlyValid,
2011                                                                  ViewCellContainer &viewCells,
2012                                                                  bool onlyUnmailed) const
2013{
2014        stack<VspNode *> nodeStack;
2015
2016        if (!root)
2017                return;
2018
2019        nodeStack.push(root);
2020       
2021        while (!nodeStack.empty())
2022        {
2023                VspNode *node = nodeStack.top();
2024                nodeStack.pop();
2025               
2026                if (node->IsLeaf())
2027                {
2028                        if (!onlyValid || node->TreeValid())
2029                        {
2030                                ViewCellLeaf *leafVc = dynamic_cast<VspLeaf *>(node)->GetViewCell();
2031
2032                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
2033                                               
2034                                if (!onlyUnmailed || !viewCell->Mailed())
2035                                {
2036                                        viewCell->Mail();
2037                                        viewCells.push_back(viewCell);
2038                                }
2039                        }
2040                }
2041                else
2042                {
2043                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2044               
2045                        nodeStack.push(interior->GetFront());
2046                        nodeStack.push(interior->GetBack());
2047                }
2048        }
2049}
2050
2051
2052int VspTree::FindNeighbors(VspLeaf *n,
2053                                                   vector<VspLeaf *> &neighbors,
2054                                                   const bool onlyUnmailed) const
2055{
2056        stack<VspNode *> nodeStack;
2057        nodeStack.push(mRoot);
2058
2059        const AxisAlignedBox3 box = GetBoundingBox(n);
2060
2061        while (!nodeStack.empty())
2062        {
2063                VspNode *node = nodeStack.top();
2064                nodeStack.pop();
2065
2066                if (node->IsLeaf())
2067                {
2068                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2069
2070                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
2071                                neighbors.push_back(leaf);
2072                }
2073                else
2074                {
2075                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2076                       
2077                        if (interior->GetPosition() > box.Max(interior->GetAxis()))
2078                                nodeStack.push(interior->GetBack());
2079                        else
2080                        {
2081                                if (interior->GetPosition() < box.Min(interior->GetAxis()))
2082                                        nodeStack.push(interior->GetFront());
2083                                else
2084                                {
2085                                        // random decision
2086                                        nodeStack.push(interior->GetBack());
2087                                        nodeStack.push(interior->GetFront());
2088                                }
2089                        }
2090                }
2091        }
2092
2093        return (int)neighbors.size();
2094}
2095
2096
2097// Find random neighbor which was not mailed
2098VspLeaf *VspTree::GetRandomLeaf(const Plane3 &plane)
2099{
2100        stack<VspNode *> nodeStack;
2101        nodeStack.push(mRoot);
2102 
2103        int mask = rand();
2104 
2105        while (!nodeStack.empty())
2106        {
2107                VspNode *node = nodeStack.top();
2108               
2109                nodeStack.pop();
2110               
2111                if (node->IsLeaf())
2112                {
2113                        return dynamic_cast<VspLeaf *>(node);
2114                }
2115                else
2116                {
2117                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2118                        VspNode *next;
2119                       
2120                        if (GetBoundingBox(interior->GetBack()).Side(plane) < 0)
2121                        {
2122                                next = interior->GetFront();
2123                        }
2124            else
2125                        {
2126                                if (GetBoundingBox(interior->GetFront()).Side(plane) < 0)
2127                                {
2128                                        next = interior->GetBack();
2129                                }
2130                                else
2131                                {
2132                                        // random decision
2133                                        if (mask & 1)
2134                                                next = interior->GetBack();
2135                                        else
2136                                                next = interior->GetFront();
2137                                        mask = mask >> 1;
2138                                }
2139                        }
2140                       
2141                        nodeStack.push(next);
2142                }
2143        }
2144 
2145        return NULL;
2146}
2147
2148
2149VspLeaf *VspTree::GetRandomLeaf(const bool onlyUnmailed)
2150{
2151        stack<VspNode *> nodeStack;
2152
2153        nodeStack.push(mRoot);
2154
2155        int mask = rand();
2156
2157        while (!nodeStack.empty())
2158        {
2159                VspNode *node = nodeStack.top();
2160                nodeStack.pop();
2161
2162                if (node->IsLeaf())
2163                {
2164                        if ( (!onlyUnmailed || !node->Mailed()) )
2165                                return dynamic_cast<VspLeaf *>(node);
2166                }
2167                else
2168                {
2169                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2170
2171                        // random decision
2172                        if (mask & 1)
2173                                nodeStack.push(interior->GetBack());
2174                        else
2175                                nodeStack.push(interior->GetFront());
2176
2177                        mask = mask >> 1;
2178                }
2179        }
2180
2181        return NULL;
2182}
2183
2184
2185void VspTree::CollectPvs(const RayInfoContainer &rays,
2186                                                 ObjectContainer &objects) const
2187{
2188        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2189
2190        Intersectable::NewMail();
2191
2192        for (rit = rays.begin(); rit != rays.end(); ++ rit)
2193        {
2194                VssRay *ray = (*rit).mRay;
2195
2196                Intersectable *object;
2197                object = ray->mOriginObject;
2198
2199        if (object)
2200                {
2201                        if (!object->Mailed())
2202                        {
2203                                object->Mail();
2204                                objects.push_back(object);
2205                        }
2206                }
2207
2208                object = ray->mTerminationObject;
2209
2210                if (object)
2211                {
2212                        if (!object->Mailed())
2213                        {
2214                                object->Mail();
2215                                objects.push_back(object);
2216                        }
2217                }
2218        }
2219}
2220
2221
2222int VspTree::ComputePvsSize(const RayInfoContainer &rays) const
2223{
2224        int pvsSize = 0;
2225
2226        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2227
2228        Intersectable::NewMail();
2229
2230        for (rit = rays.begin(); rit != rays.end(); ++ rit)
2231        {
2232                VssRay *ray = (*rit).mRay;
2233
2234                if (ray->mOriginObject)
2235                {
2236                        if (!ray->mOriginObject->Mailed())
2237                        {
2238                                ray->mOriginObject->Mail();
2239                                ++ pvsSize;
2240                        }
2241                }
2242                if (ray->mTerminationObject)
2243                {
2244                        if (!ray->mTerminationObject->Mailed())
2245                        {
2246                                ray->mTerminationObject->Mail();
2247                                ++ pvsSize;
2248                        }
2249                }
2250        }
2251
2252        return pvsSize;
2253}
2254
2255
2256float VspTree::GetEpsilon() const
2257{
2258        return mEpsilon;
2259}
2260
2261
2262int VspTree::CastLineSegment(const Vector3 &origin,
2263                                                         const Vector3 &termination,
2264                             ViewCellContainer &viewcells)
2265{
2266        int hits = 0;
2267
2268        float mint = 0.0f, maxt = 1.0f;
2269        const Vector3 dir = termination - origin;
2270
2271        stack<LineTraversalData> tStack;
2272
2273        Intersectable::NewMail();
2274        //ViewCell::NewMail();
2275
2276        Vector3 entp = origin;
2277        Vector3 extp = termination;
2278
2279        VspNode *node = mRoot;
2280        VspNode *farChild;
2281
2282        float position;
2283        int axis;
2284
2285        while (1)
2286        {
2287                if (!node->IsLeaf())
2288                {
2289                        VspInterior *in = dynamic_cast<VspInterior *>(node);
2290                        position = in->GetPosition();
2291                        axis = in->GetAxis();
2292
2293                        if (entp[axis] <= position)
2294                        {
2295                                if (extp[axis] <= position)
2296                                {
2297                                        node = in->GetBack();
2298                                        // cases N1,N2,N3,P5,Z2,Z3
2299                                        continue;
2300                                } else
2301                                {
2302                                        // case N4
2303                                        node = in->GetBack();
2304                                        farChild = in->GetFront();
2305                                }
2306                        }
2307                        else
2308                        {
2309                                if (position <= extp[axis])
2310                                {
2311                                        node = in->GetFront();
2312                                        // cases P1,P2,P3,N5,Z1
2313                                        continue;
2314                                }
2315                                else
2316                                {
2317                                        node = in->GetFront();
2318                                        farChild = in->GetBack();
2319                                        // case P4
2320                                }
2321                        }
2322
2323                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2324                        // case N4 or P4
2325                        const float tdist = (position - origin[axis]) / dir[axis];
2326                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2327
2328                        extp = origin + dir * tdist;
2329                        maxt = tdist;
2330                }
2331                else
2332                {
2333                        // compute intersection with all objects in this leaf
2334                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2335                        ViewCell *vc = leaf->GetViewCell();
2336
2337                        // don't have to mail because each view cell belongs to exactly one leaf
2338                        //if (!vc->Mailed())
2339                        //{
2340                        //      vc->Mail();
2341                                viewcells.push_back(vc);
2342                                ++ hits;
2343                        //}
2344#if 0
2345                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2346#endif
2347                        // get the next node from the stack
2348                        if (tStack.empty())
2349                                break;
2350
2351                        entp = extp;
2352                        mint = maxt;
2353                       
2354                        LineTraversalData &s  = tStack.top();
2355                        node = s.mNode;
2356                        extp = s.mExitPoint;
2357                        maxt = s.mMaxT;
2358                        tStack.pop();
2359                }
2360        }
2361
2362        return hits;
2363}
2364
2365
2366int VspTree::CastRay(Ray &ray)
2367{
2368        int hits = 0;
2369
2370        stack<LineTraversalData> tStack;
2371        const Vector3 dir = ray.GetDir();
2372
2373        float maxt, mint;
2374
2375        if (!mBoundingBox.GetRaySegment(ray, mint, maxt))
2376                return 0;
2377
2378        Intersectable::NewMail();
2379        ViewCell::NewMail();
2380
2381        Vector3 entp = ray.Extrap(mint);
2382        Vector3 extp = ray.Extrap(maxt);
2383
2384        const Vector3 origin = entp;
2385
2386        VspNode *node = mRoot;
2387        VspNode *farChild = NULL;
2388
2389        float position;
2390        int axis;
2391
2392        while (1)
2393        {
2394                if (!node->IsLeaf())
2395                {
2396                        VspInterior *in = dynamic_cast<VspInterior *>(node);
2397                        position = in->GetPosition();
2398                        axis = in->GetAxis();
2399
2400                        if (entp[axis] <= position)
2401                        {
2402                                if (extp[axis] <= position)
2403                                {
2404                                        node = in->GetBack();
2405                                        // cases N1,N2,N3,P5,Z2,Z3
2406                                        continue;
2407                                }
2408                                else
2409                                {
2410                                        // case N4
2411                                        node = in->GetBack();
2412                                        farChild = in->GetFront();
2413                                }
2414                        }
2415                        else
2416                        {
2417                                if (position <= extp[axis])
2418                                {
2419                                        node = in->GetFront();
2420                                        // cases P1,P2,P3,N5,Z1
2421                                        continue;
2422                                }
2423                                else
2424                                {
2425                                        node = in->GetFront();
2426                                        farChild = in->GetBack();
2427                                        // case P4
2428                                }
2429                        }
2430
2431                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2432                        // case N4 or P4
2433                        const float tdist = (position - origin[axis]) / dir[axis];
2434                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2435                        extp = origin + dir * tdist;
2436                        maxt = tdist;
2437                }
2438                else
2439                {
2440                        // compute intersection with all objects in this leaf
2441                        VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2442                        ViewCell *vc = leaf->GetViewCell();
2443
2444                        if (!vc->Mailed())
2445                        {
2446                                vc->Mail();
2447                                // todo: add view cells to ray
2448                                ++ hits;
2449                        }
2450#if 0
2451                        leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2452#endif
2453                        // get the next node from the stack
2454                        if (tStack.empty())
2455                                break;
2456
2457                        entp = extp;
2458                        mint = maxt;
2459                       
2460                        LineTraversalData &s  = tStack.top();
2461                        node = s.mNode;
2462                        extp = s.mExitPoint;
2463                        maxt = s.mMaxT;
2464                        tStack.pop();
2465                }
2466        }
2467
2468        return hits;
2469}
2470
2471
2472ViewCell *VspTree::GetViewCell(const Vector3 &point, const bool active)
2473{
2474        if (mRoot == NULL)
2475                return NULL;
2476
2477        stack<VspNode *> nodeStack;
2478        nodeStack.push(mRoot);
2479 
2480        ViewCellLeaf *viewcell = NULL;
2481 
2482        while (!nodeStack.empty()) 
2483        {
2484                VspNode *node = nodeStack.top();
2485                nodeStack.pop();
2486       
2487                if (node->IsLeaf())
2488                {
2489                        viewcell = dynamic_cast<VspLeaf *>(node)->GetViewCell();
2490                        break;
2491                }
2492                else   
2493                {       
2494                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2495     
2496                        // random decision
2497                        if (interior->GetPosition() - point[interior->GetAxis()] < 0)
2498                                nodeStack.push(interior->GetBack());
2499                        else
2500                                nodeStack.push(interior->GetFront());
2501                }
2502        }
2503 
2504        if (active)
2505                return mViewCellsTree->GetActiveViewCell(viewcell);
2506        else
2507                return viewcell;
2508}
2509
2510
2511bool VspTree::ViewPointValid(const Vector3 &viewPoint) const
2512{
2513        VspNode *node = mRoot;
2514
2515        while (1)
2516        {
2517                // early exit
2518                if (node->TreeValid())
2519                        return true;
2520
2521                if (node->IsLeaf())
2522                        return false;
2523                       
2524                VspInterior *in = dynamic_cast<VspInterior *>(node);
2525                                       
2526                if (in->GetPosition() - viewPoint[in->GetAxis()] <= 0)
2527                {
2528                        node = in->GetBack();
2529                }
2530                else
2531                {
2532                        node = in->GetFront();
2533                }
2534        }
2535
2536        // should never come here
2537        return false;
2538}
2539
2540
2541void VspTree::PropagateUpValidity(VspNode *node)
2542{
2543        const bool isValid = node->TreeValid();
2544
2545        // propagative up invalid flag until only invalid nodes exist over this node
2546        if (!isValid)
2547        {
2548                while (!node->IsRoot() && node->GetParent()->TreeValid())
2549                {
2550                        node = node->GetParent();
2551                        node->SetTreeValid(false);
2552                }
2553        }
2554        else
2555        {
2556                // propagative up valid flag until one of the subtrees is invalid
2557                while (!node->IsRoot() && !node->TreeValid())
2558                {
2559            node = node->GetParent();
2560                        VspInterior *interior = dynamic_cast<VspInterior *>(node);
2561                       
2562                        // the parent is valid iff both leaves are valid
2563                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
2564                                                           interior->GetFront()->TreeValid());
2565                }
2566        }
2567}
2568
2569#if ZIPPED_VIEWCELLS
2570bool VspTree::Export(ogzstream &stream)
2571#else
2572bool VspTree::Export(ofstream &stream)
2573#endif
2574{
2575        ExportNode(mRoot, stream);
2576
2577        return true;
2578}
2579
2580
2581#if ZIPPED_VIEWCELLS
2582void VspTree::ExportNode(VspNode *node, ogzstream &stream)
2583#else
2584void VspTree::ExportNode(VspNode *node, ofstream &stream)
2585#endif
2586{
2587        if (node->IsLeaf())
2588        {
2589                VspLeaf *leaf = dynamic_cast<VspLeaf *>(node);
2590                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
2591
2592                int id = -1;
2593                if (viewCell != mOutOfBoundsCell)
2594                        id = viewCell->GetId();
2595
2596                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
2597        }
2598        else
2599        {       
2600                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2601       
2602                AxisAlignedPlane plane = interior->GetPlane();
2603                stream << "<Interior plane=\"" << plane.mPosition << " "
2604                           << plane.mAxis << "\">" << endl;
2605
2606                ExportNode(interior->GetBack(), stream);
2607                ExportNode(interior->GetFront(), stream);
2608
2609                stream << "</Interior>" << endl;
2610        }
2611}
2612
2613
2614int VspTree::SplitRays(const AxisAlignedPlane &plane,
2615                                           RayInfoContainer &rays,
2616                                           RayInfoContainer &frontRays,
2617                                           RayInfoContainer &backRays) const
2618{
2619        int splits = 0;
2620
2621        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2622
2623        for (rit = rays.begin(); rit != rit_end; ++ rit)
2624        {
2625                RayInfo bRay = *rit;
2626               
2627                VssRay *ray = bRay.mRay;
2628                float t;
2629
2630                // get classification and receive new t
2631                //-- test if start point behind or in front of plane
2632                const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
2633                       
2634                if (side == 0)
2635                {
2636                        ++ splits;
2637
2638                        if (ray->HasPosDir(plane.mAxis))
2639                        {
2640                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2641                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2642                        }
2643                        else
2644                        {
2645                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2646                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2647                        }
2648                }
2649                else if (side == 1)
2650                {
2651                        frontRays.push_back(bRay);
2652                }
2653                else
2654                {
2655                        backRays.push_back(bRay);
2656                }
2657        }
2658
2659        return splits;
2660}
2661
2662
2663AxisAlignedBox3 VspTree::GetBoundingBox(VspNode *node) const
2664{
2665        if (!node->GetParent())
2666                return mBoundingBox;
2667
2668        if (!node->IsLeaf())
2669        {
2670                return (dynamic_cast<VspInterior *>(node))->GetBoundingBox();           
2671        }
2672
2673        VspInterior *parent = dynamic_cast<VspInterior *>(node->GetParent());
2674
2675        AxisAlignedBox3 box(parent->GetBoundingBox());
2676
2677        if (parent->GetFront() == node)
2678      box.SetMin(parent->GetAxis(), parent->GetPosition());
2679    else
2680      box.SetMax(parent->GetAxis(), parent->GetPosition());
2681
2682        return box;
2683}
2684
2685
2686int VspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box,
2687                                                                         ViewCellContainer &viewCells) const
2688{
2689        stack<VspNode *> nodeStack;
2690 
2691        ViewCell::NewMail();
2692
2693        while (!nodeStack.empty())
2694        {
2695                VspNode *node = nodeStack.top();
2696                nodeStack.pop();
2697
2698                const AxisAlignedBox3 bbox = GetBoundingBox(node);
2699
2700                if (bbox.Includes(box))
2701                {
2702                        // node geometry is contained in box
2703                        CollectViewCells(node, true, viewCells, true);
2704                }
2705                else if (Overlap(bbox, box))
2706                {
2707                        if (node->IsLeaf())
2708                        {
2709                                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2710                       
2711                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
2712                                {
2713                                        leaf->GetViewCell()->Mail();
2714                                        viewCells.push_back(leaf->GetViewCell());
2715                                }
2716                        }
2717                        else
2718                        {
2719                                VspInterior *interior = dynamic_cast<VspInterior *>(node);
2720                       
2721                                VspNode *first = interior->GetFront();
2722                                VspNode *second = interior->GetBack();
2723           
2724                                nodeStack.push(first);
2725                                nodeStack.push(second);
2726                        }
2727                }       
2728                // default: cull
2729        }
2730
2731        return (int)viewCells.size();
2732}
2733
2734
2735void VspTree::CollectDirtyCandidates(VspSplitCandidate *sc,
2736                                                                         vector<SplitCandidate *> &dirtyList)
2737{
2738        VspTraversalData &tData = sc->mParentData;
2739        VspLeaf *node = tData.mNode;
2740       
2741        KdLeaf::NewMail();
2742
2743        RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
2744
2745        // add all kd nodes seen by the rays
2746        for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
2747        {
2748                VssRay *ray = (*rit).mRay;
2749       
2750                KdLeaf *leaf = mOspTree->GetLeaf(ray->mOrigin, ray->mOriginNode);
2751       
2752                if (!leaf->Mailed())
2753                {
2754                        leaf->Mail();
2755                        dirtyList.push_back(leaf->mSplitCandidate);
2756                }
2757               
2758                leaf = mOspTree->GetLeaf(ray->mTermination, ray->mTerminationNode);
2759       
2760                if (!leaf->Mailed())
2761                {
2762                        leaf->Mail();
2763                        dirtyList.push_back(leaf->mSplitCandidate);
2764                }
2765        }
2766}
2767
2768
2769void VspTree::ProcessRays(const VssRayContainer &sampleRays,
2770                                                  RayInfoContainer &rays)
2771{
2772        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
2773
2774        long startTime = GetTime();
2775
2776        cout << "storing rays ... ";
2777
2778        Intersectable::NewMail();
2779
2780        //-- store rays and objects
2781        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
2782        {
2783                VssRay *ray = *rit;
2784                float minT, maxT;
2785                static Ray hray;
2786
2787                hray.Init(*ray);
2788               
2789                // TODO: not very efficient to implictly cast between rays types
2790                if (GetBoundingBox().GetRaySegment(hray, minT, maxT))
2791                {
2792                        float len = ray->Length();
2793
2794                        if (!len)
2795                                len = Limits::Small;
2796
2797                        rays.push_back(RayInfo(ray, minT / len, maxT / len));
2798                }
2799        }
2800
2801        cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
2802}
2803
2804
2805void VspTree::GetViewCells(const VssRay &ray, ViewCellContainer &viewCells)
2806{
2807        // if no precomputation of view cells
2808        CastLineSegment(ray.mOrigin, ray.mTermination, viewCells);
2809}
2810
2811
2812/*******************************************************************/
2813/*                  class OspTree implementation                   */
2814/*******************************************************************/
2815
2816
2817OspTree::OspTree():
2818mRoot(NULL),
2819mTimeStamp(1),
2820mCopyFromKdTree(false)
2821{
2822        ReadEnvironment();
2823
2824        mSplitCandidates = new vector<SortableEntry>;
2825}
2826
2827
2828OspTree::OspTree(const KdTree &kdTree)
2829{
2830        ReadEnvironment();
2831
2832        // copy values from kd tree
2833        mCopyFromKdTree = true;
2834        mBoundingBox = kdTree.GetBox();
2835        mRoot = kdTree.GetRoot();
2836}
2837
2838
2839OspTree::~OspTree()
2840{
2841        KdIntersectableMap::iterator it, it_end = mKdIntersectables.end();
2842
2843        for (it = mKdIntersectables.begin(); it != mKdIntersectables.end(); ++ it)
2844        {
2845                DEL_PTR((*it).second);
2846        }
2847
2848        // if not using kd tree root
2849        if (!mCopyFromKdTree)
2850                DEL_PTR(mRoot);
2851}
2852
2853
2854void OspTree::ReadEnvironment()
2855{
2856        bool randomize = false;
2857        Environment::GetSingleton()->GetBoolValue("VspTree.Construction.randomize", randomize);
2858        if (randomize)
2859                Randomize(); // initialise random generator for heuristics
2860
2861        //-- termination criteria for autopartition
2862        Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxDepth", mTermMaxDepth);
2863        Environment::GetSingleton()->GetIntValue("OspTree.Termination.maxLeaves", mTermMaxLeaves);
2864        Environment::GetSingleton()->GetIntValue("OspTree.Termination.minObjects", mTermMinObjects);
2865        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minProbability", mTermMinProbability);
2866       
2867        Environment::GetSingleton()->GetIntValue("OspTree.Termination.missTolerance", mTermMissTolerance);
2868       
2869        //-- max cost ratio for early tree termination
2870        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.maxCostRatio", mTermMaxCostRatio);
2871
2872        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.minGlobalCostRatio",
2873                mTermMinGlobalCostRatio);
2874       
2875
2876        //-- factors for bsp tree split plane heuristics
2877        Environment::GetSingleton()->GetFloatValue("OspTree.Termination.ct_div_ci", mCtDivCi);
2878
2879        //-- partition criteria
2880        Environment::GetSingleton()->GetFloatValue("OspTree.Construction.epsilon", mEpsilon);
2881
2882        // if only the driving axis is used for axis aligned split
2883        Environment::GetSingleton()->GetBoolValue("OspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
2884        Environment::GetSingleton()->GetFloatValue("OspTree.maxStaticMemory", mMaxMemory);
2885        Environment::GetSingleton()->GetBoolValue("OspTree.useCostHeuristics", mUseCostHeuristics);
2886
2887
2888        char subdivisionStatsLog[100];
2889        Environment::GetSingleton()->GetStringValue("OspTree.subdivisionStats", subdivisionStatsLog);
2890        mSubdivisionStats.open(subdivisionStatsLog);
2891
2892        Environment::GetSingleton()->GetFloatValue("OspTree.Construction.splitBorder", mSplitBorder);
2893        Environment::GetSingleton()->GetFloatValue("OspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
2894       
2895
2896        //-- debug output
2897
2898        Debug << "******* OSP tree options ******** " << endl;
2899
2900    Debug << "max depth: " << mTermMaxDepth << endl;
2901        Debug << "min probabiliy: " << mTermMinProbability << endl;
2902        Debug << "min objects: " << mTermMinObjects << endl;
2903        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
2904        Debug << "miss tolerance: " << mTermMissTolerance << endl;
2905        Debug << "max leaves: " << mTermMaxLeaves << endl;
2906
2907        Debug << "randomize: " << randomize << endl;
2908
2909        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
2910        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
2911        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
2912        Debug << "max memory: " << mMaxMemory << endl;
2913        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
2914        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
2915       
2916        Debug << "split borders: " << mSplitBorder << endl;
2917        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
2918        Debug << endl;
2919}
2920
2921
2922void OspTree::SplitObjects(KdLeaf *parent,
2923                                                   const AxisAlignedPlane& splitPlane,
2924                                                   const ObjectContainer &objects,
2925                                                   ObjectContainer &front,
2926                                                   ObjectContainer &back)
2927{
2928        ObjectContainer::const_iterator oit, oit_end = objects.end();
2929
2930    for (oit = objects.begin(); oit != oit_end; ++ oit)
2931        {
2932                Intersectable *object = *oit;
2933               
2934                // initialise
2935                if (parent->IsRoot())
2936                        object->mReferences = 1;
2937
2938                // remove parent ref
2939                -- object->mReferences;
2940
2941                // determine the side of this ray with respect to the plane
2942                const AxisAlignedBox3 box = object->GetBox();
2943
2944                if (box.Max(splitPlane.mAxis) >= splitPlane.mPosition)
2945                {
2946            front.push_back(object);
2947                        ++ object->mReferences;
2948                }
2949
2950                if (box.Min(splitPlane.mAxis) < splitPlane.mPosition)
2951                {
2952                        back.push_back(object);
2953                        ++ object->mReferences;
2954                }
2955        }
2956
2957        mOspStats.objectRefs -= (int)objects.size();
2958        mOspStats.objectRefs += (int)back.size() + (int)front.size();
2959}
2960
2961
2962KdInterior *OspTree::SubdivideNode(
2963                                                                   const AxisAlignedPlane &splitPlane,
2964                                                                   const OspTraversalData &tData,
2965                                                                   OspTraversalData &frontData,
2966                                                                   OspTraversalData &backData
2967                                                                   )
2968{
2969        KdLeaf *leaf = tData.mNode;
2970       
2971        // two new leaves
2972    mOspStats.nodes += 2;
2973
2974        // add the new nodes to the tree
2975        KdInterior *node = new KdInterior(leaf->mParent);
2976
2977        const int axis = splitPlane.mAxis;
2978        const float position = splitPlane.mPosition;
2979
2980        node->mAxis = axis;
2981        node->mPosition = position;
2982        node->mBox = tData.mBoundingBox;
2983
2984
2985        //-- the front and back traversal data is filled with the new values
2986
2987        // bounding boxes: split front and back node geometry
2988        tData.mBoundingBox.Split(splitPlane.mAxis, splitPlane.mPosition,
2989                frontData.mBoundingBox, backData.mBoundingBox);
2990
2991        frontData.mDepth = backData.mDepth = tData.mDepth + 1;
2992               
2993        // TODO matt: compute pvs
2994        frontData.mPvs = 999; // ComputePvsSize(*frontData.mRays);
2995        backData.mPvs = 999; //ComputePvsSize(*backData.mRays);
2996
2997        frontData.mProbability = frontData.mBoundingBox.GetVolume();
2998
2999       
3000        //-- subdivide rays
3001        frontData.mRays = new RayInfoContainer();
3002        backData.mRays = new RayInfoContainer();
3003
3004        SplitRays(splitPlane,
3005                          *tData.mRays,
3006                          *frontData.mRays,
3007                          *backData.mRays);
3008
3009
3010        //-- classify objects
3011        int objectsBack = 0;
3012        int objectsFront = 0;
3013
3014        ObjectContainer::const_iterator mi, mi_end = leaf->mObjects.end();
3015
3016        for ( mi = leaf->mObjects.begin(); mi != mi_end; ++ mi)
3017        {
3018                // determine the side of this ray with respect to the plane
3019                const AxisAlignedBox3 box = (*mi)->GetBox();
3020               
3021                if (box.Max(axis) > position + mEpsilon)
3022                        ++ objectsFront;
3023   
3024                if (box.Min(axis) < position - mEpsilon)
3025                        ++ objectsBack;
3026        }
3027
3028        KdLeaf *back = new KdLeaf(node, objectsBack);
3029        KdLeaf *front = new KdLeaf(node, objectsFront);
3030
3031        /////////////
3032        //-- create front and back leaf
3033
3034        KdInterior *parent = leaf->mParent;
3035
3036        // replace a link from node's parent
3037        if (parent)
3038        {
3039                parent->ReplaceChildLink(leaf, node);
3040                node->mParent = parent;
3041        }
3042        else // new root
3043        {
3044                mRoot = node;
3045        }
3046
3047        // and setup child links
3048        node->SetupChildLinks(back, front);
3049
3050        SplitObjects(leaf, splitPlane, leaf->mObjects, front->mObjects, back->mObjects);
3051
3052        ProcessMultipleRefs(front);
3053    ProcessMultipleRefs(back);
3054
3055        backData.mNode = back;
3056        frontData.mNode = front;
3057
3058        //delete leaf;
3059        return node;
3060}
3061
3062
3063KdNode *OspTree::Subdivide(SplitQueue &tQueue,
3064                                                   SplitCandidate *splitCandidate,
3065                                                   const bool globalCriteriaMet)
3066{
3067        OspSplitCandidate *sc = dynamic_cast<OspSplitCandidate *>(splitCandidate);
3068        OspTraversalData &tData = sc->mParentData;
3069
3070        KdNode *newNode = tData.mNode;
3071
3072        if (!LocalTerminationCriteriaMet(tData) && !globalCriteriaMet)
3073        {       
3074                OspTraversalData tFrontData;
3075                OspTraversalData tBackData;
3076
3077                //-- continue subdivision
3078               
3079                // create new interior node and two leaf node
3080                const AxisAlignedPlane splitPlane = sc->mSplitPlane;
3081               
3082                newNode = SubdivideNode(splitPlane,
3083                                                                tData,
3084                                                                tFrontData,
3085                                                                tBackData);
3086       
3087                const int maxCostMisses = sc->mMaxCostMisses;
3088
3089                // how often was max cost ratio missed in this branch?
3090                tFrontData.mMaxCostMisses = maxCostMisses;
3091                tBackData.mMaxCostMisses = maxCostMisses;
3092                       
3093                //-- push the new split candidates on the queue
3094                OspSplitCandidate *frontCandidate = new OspSplitCandidate(tFrontData);
3095                OspSplitCandidate *backCandidate = new OspSplitCandidate(tBackData);
3096
3097                EvalSplitCandidate(*frontCandidate);
3098                EvalSplitCandidate(*backCandidate);
3099
3100                tQueue.Push(frontCandidate);
3101                tQueue.Push(backCandidate);
3102
3103                // delete old leaf node
3104                DEL_PTR(tData.mNode);
3105        }
3106
3107
3108        //-- terminate traversal
3109        if (newNode->IsLeaf())
3110        {
3111                EvaluateLeafStats(tData);               
3112        }
3113       
3114        //-- cleanup
3115        tData.Clear();
3116       
3117        return newNode;
3118}
3119
3120
3121void OspTree::EvalSplitCandidate(OspSplitCandidate &splitCandidate)
3122{
3123        float frontProb;
3124        float backProb;
3125       
3126        // compute locally best split plane
3127        const bool success =
3128                SelectSplitPlane(splitCandidate.mParentData, splitCandidate.mSplitPlane, frontProb, backProb);
3129
3130        float oldRenderCost;
3131
3132        // compute global decrease in render cost
3133        const float renderCostDecr = EvalRenderCostDecrease(splitCandidate.mSplitPlane,
3134                                                                                                                splitCandidate.mParentData,
3135                                                                                                                oldRenderCost);
3136
3137        splitCandidate.SetRenderCostDecrease(renderCostDecr);
3138
3139#if 0
3140        const float priority = (float)-data.mDepth;
3141#else   
3142        // take render cost of node into account
3143        // otherwise danger of being stuck in a local minimum!!
3144        const float factor = mRenderCostDecreaseWeight;
3145        const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;
3146
3147#endif
3148
3149        // compute global decrease in render cost
3150        splitCandidate.SetPriority(priority);
3151
3152        splitCandidate.mMaxCostMisses =
3153                success ? splitCandidate.mParentData.mMaxCostMisses : splitCandidate.mParentData.mMaxCostMisses + 1;
3154}
3155
3156
3157bool OspTree::LocalTerminationCriteriaMet(const OspTraversalData &data) const
3158{
3159        // matt: TODO
3160        return (
3161                //(data.mNode->mObjects.size() < mTermMinObjects) ||
3162                //(data.mProbability <= mTermMinProbability) ||
3163                (data.mDepth >= mTermMaxDepth)
3164                 );
3165}
3166
3167
3168bool OspTree::GlobalTerminationCriteriaMet(const OspTraversalData &data) const
3169{
3170        // matt: TODO
3171        return (
3172                (mOspStats.Leaves() >= mTermMaxLeaves)
3173                //mOutOfMemory ||
3174                //(mGlobalCostMisses >= mTermGlobalCostMissTolerance)
3175                );
3176}
3177
3178
3179void OspTree::EvaluateLeafStats(const OspTraversalData &data)
3180{
3181        // the node became a leaf -> evaluate stats for leafs
3182        KdLeaf *leaf = data.mNode;
3183
3184        if (data.mPvs > mOspStats.maxPvs)
3185        {
3186                mOspStats.maxPvs = data.mPvs;
3187        }
3188
3189        mOspStats.pvs += data.mPvs;
3190
3191        if (data.mDepth < mOspStats.minDepth)
3192        {
3193                mOspStats.minDepth = data.mDepth;
3194        }
3195       
3196        if (data.mDepth >= mTermMaxDepth)
3197        {
3198        ++ mOspStats.maxDepthNodes;
3199                //Debug << "new max depth: " << mVspStats.maxDepthNodes << endl;
3200        }
3201
3202//      if (data.mPvs < mTermMinPvs)
3203//              ++ mOspStats.minPvsNodes;
3204
3205        if (data.mProbability <= mTermMinProbability)
3206                ++ mOspStats.minProbabilityNodes;
3207       
3208        // accumulate depth to compute average depth
3209        mOspStats.accumDepth += data.mDepth;
3210
3211        ++ mCreatedLeaves;
3212
3213#ifdef _DEBUG
3214        Debug << "BSP stats: "
3215                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
3216                 // << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
3217                  << "Prob: " << data.mProbability << " (min: " << mTermMinProbability << "), "
3218                  << "#pvs: " << data.mPvs << ")\n";
3219#endif
3220}
3221
3222
3223float OspTree::EvalLocalCostHeuristics(const OspTraversalData &tData,
3224                                                                           const AxisAlignedBox3 &box,
3225                                                                           const int axis,
3226                                                                           float &position,
3227                                                                           int &objectsFront,
3228                                                                           int &objectsBack)
3229{
3230        // go through the lists, count the number of objects left and right
3231        // and evaluate the following cost funcion:
3232        // C = ct_div_ci  + (ol + or)/queries
3233       
3234        const float minBox = box.Min(axis);
3235        const float maxBox = box.Max(axis);
3236
3237        const float sizeBox = maxBox - minBox;
3238
3239        float minBand = minBox + mSplitBorder * (maxBox - minBox);
3240        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
3241
3242        //-- sort so we can use a sweep
3243        SortSplitCandidates(tData, axis, minBand, maxBand);
3244
3245        int numViewCells;
3246
3247        float totalVol = PrepareHeuristics(tData, numViewCells);
3248        float voll = 0;
3249        float volr = totalVol;
3250
3251        int vcl = 0;
3252        int vcr = numViewCells;
3253 
3254        const int totalPvs = tData.mNode->mObjects.size();
3255       
3256        int pvsl = 0;
3257        int pvsr = totalPvs;
3258
3259        float sum = (float)totalVol * sizeBox;
3260
3261        /////////////////////////////////
3262
3263        // note: initialised to take mid split
3264        // if no good split can be found,
3265        position = minBox + 0.5f * sizeBox;
3266       
3267        // the relative cost ratio
3268        float ratio = 99999999.0f;
3269        bool splitPlaneFound = false;
3270
3271        float volBack = voll;
3272        float volFront = volr;
3273
3274        int pvsBack = pvsl;
3275        int pvsFront = pvsr;
3276       
3277        float minSum = 1e20f;
3278
3279        /////////////////////////////
3280
3281        Intersectable::NewMail();
3282        ViewCell::NewMail();
3283
3284        //-- traverse through events and find best split plane
3285
3286        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
3287
3288        for (ci = mSplitCandidates->begin(); ci != ci_end; ++ ci)
3289        {
3290                Intersectable *object = (*ci).mObject;
3291
3292                EvalHeuristicsContribution(tData.mNode,
3293                                                                   *ci,
3294                                                                   voll, volr,
3295                                                                   pvsl, pvsr,
3296                                                                   vcl, vcr);
3297
3298                // Note: sufficient to compare size of bounding boxes of front and back side?
3299
3300                if (((*ci).mPos >= minBand) && ((*ci).mPos <= maxBand))
3301                {
3302                        sum = voll * pvsl + volr * pvsr;
3303
3304                        float currentPos;
3305                       
3306                        // HACK: current positition is BETWEEN visibility events
3307                        if (0 && ((ci + 1) != ci_end))
3308                                currentPos = ((*ci).mPos + (*(ci + 1)).mPos) * 0.5f;
3309                        else
3310                                currentPos = (*ci).mPos;                       
3311
3312                        // note matt: can happen that volr is less than zero: bug or numerical error?
3313                        //if (volr < 0) Debug << "warning!";
3314
3315                        /*Debug << "pos: " << (*ci).mPos
3316                                 << "\t (pvsl: " << pvsl << ", pvsr: " << pvsr << ")"
3317                                 << "\t (voll: " << voll << ", volr: " << volr << ")"
3318                                 << "\t (vcl: " << vcl << ", vcr: " << vcr << ", nvc: " << numViewCells << ")"
3319                                 << "\t sum: " << sum << endl;*/
3320
3321                        if (sum < minSum)
3322                        {
3323                                splitPlaneFound = true;
3324
3325                                minSum = sum;
3326
3327                                pvsBack = pvsl;
3328                                pvsFront = pvsr;
3329                               
3330                                volBack = voll;
3331                                volFront = volr;
3332
3333                                position = currentPos;
3334                        }
3335                }
3336        }
3337       
3338
3339        //-- compute cost
3340
3341        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
3342        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
3343
3344        const float pOverall = totalVol;
3345        const float pBack = volBack;
3346        const float pFront = volFront;
3347
3348        const float penaltyOld = (float)totalPvs;
3349    const float penaltyFront = (float)pvsFront;
3350        const float penaltyBack = (float)pvsBack;
3351       
3352        const float oldRenderCost = penaltyOld * pOverall + Limits::Small;
3353        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
3354
3355
3356        if (splitPlaneFound)
3357        {
3358                ratio = newRenderCost / oldRenderCost;
3359        }
3360       
3361        //Debug << "axis=" << axis << " costRatio=" << ratio << " pos="
3362        //        << position << " t=" << (position - minBox) / (maxBox - minBox)
3363        //      << "\t pb=(" << volBack << ")\t pf=(" << volFront << ")" << endl;
3364
3365        const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();
3366
3367        Debug << "\n§§§§ eval local cost §§§§" << endl
3368                  << "back pvs: " << pvsBack << " front pvs: " << pvsFront << " total pvs: " << totalPvs << endl
3369                  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol << " p: " << pOverall / viewSpaceVol << endl
3370                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
3371                  << "render cost decrease: " << oldRenderCost / viewSpaceVol - newRenderCost / viewSpaceVol << endl;
3372
3373        return ratio;
3374}
3375
3376
3377void OspTree::SortSplitCandidates(const OspTraversalData &tData,
3378                                                                  const int axis,
3379                                                                  float minBand,
3380                                                                  float maxBand)
3381
3382{
3383        mSplitCandidates->clear();
3384
3385        RayInfoContainer *rays = tData.mRays;
3386        KdLeaf *leaf = tData.mNode;
3387
3388        int requestedSize = 2 * (int)rays->size() + 2 * (int)leaf->mObjects.size();
3389
3390        // creates a sorted split candidates array
3391        if (mSplitCandidates->capacity() > 500000 &&
3392                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
3393        {
3394        delete mSplitCandidates;
3395                mSplitCandidates = new vector<SortableEntry>;
3396        }
3397
3398        mSplitCandidates->reserve(requestedSize);
3399
3400        float pos;
3401
3402        //-- insert ray queries
3403        //-- we are intersested only in rays which intersect an object that
3404        //-- is part of the kd node because they can induce a change in view cell
3405        //-- volume on the left and the right part.
3406        //-- Note that they cannot induce a change in pvs size!!
3407
3408        RayInfoContainer::const_iterator rit, rit_end = rays->end();
3409
3410        for (rit = rays->begin(); rit < rit_end; ++ rit)
3411        {
3412                VssRay *ray = (*rit).mRay;
3413
3414                const bool positive = ray->HasPosDir(axis);
3415               
3416                // if hitpoint with object is inside this node
3417                if (ray->mOriginObject && (GetLeaf(ray->mOrigin, ray->mOriginNode) == leaf))
3418                {
3419                        pos = ray->mOrigin[axis];
3420       
3421                        mSplitCandidates->push_back(
3422                                SortableEntry(SortableEntry::BOX_INTERSECT,
3423                                                          pos,
3424                                                          ray->mOriginObject,
3425                                                          ray)
3426                                                          );
3427                }
3428
3429                if (ray->mTerminationObject && (GetLeaf(ray->mTermination, ray->mTerminationNode) == leaf))
3430                {
3431                        pos = ray->mTermination[axis];
3432
3433                        mSplitCandidates->push_back(
3434                                SortableEntry(SortableEntry::BOX_INTERSECT,
3435                                                          pos,
3436                                                          ray->mTerminationObject,
3437                                                          ray)
3438                                                          );
3439                }
3440        }
3441
3442
3443        //-- insert object queries
3444        //-- These queries can induce a change in pvs size.
3445        //-- Note that they cannot induce a change in view cell volume, as
3446        //-- there is no ray associated with these events!!
3447
3448        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
3449
3450    for (oit = leaf->mObjects.begin(); oit != leaf->mObjects.end(); ++ oit )
3451        {
3452                Intersectable *object = *oit;
3453                const AxisAlignedBox3 box = object->GetBox();
3454
3455                mSplitCandidates->push_back(
3456                        SortableEntry(SortableEntry::BOX_MIN,
3457                                                  box.Min(axis),
3458                                                  object,
3459                                                  NULL)
3460                                                  );
3461   
3462   
3463                mSplitCandidates->push_back(
3464                        SortableEntry(SortableEntry::BOX_MAX,
3465                                                  box.Max(axis),
3466                                                  object,
3467                                                  NULL)
3468                                                  );
3469        }
3470
3471        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
3472}
3473
3474
3475const OspTreeStatistics &OspTree::GetStatistics() const
3476{
3477        return mOspStats;
3478}
3479
3480
3481float OspTree::PrepareHeuristics(const VssRay &ray, int &numViewCells)
3482{
3483        float vol = 0;
3484        numViewCells = 0;
3485
3486        ViewCellContainer viewCells;
3487
3488        mVspTree->GetViewCells(ray, viewCells);
3489       
3490        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
3491
3492        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
3493        {
3494                ViewCell *vc = (*vit);
3495
3496                if (!vc->Mailed())
3497                {
3498                        //Debug << "single vol: "<< vc->GetVolume() << endl;
3499                        vc->Mail();
3500                        vc->mCounter = 0;
3501                        vol += vc->GetVolume();
3502                        ++ numViewCells;
3503                }
3504               
3505                // view cell volume already added => just increase counter
3506                ++ vc->mCounter;
3507        }
3508
3509        return vol;
3510}
3511
3512
3513float OspTree::PrepareHeuristics(const OspTraversalData &tData, int &numViewCells)
3514{       
3515        float vol = 0;
3516
3517        Intersectable::NewMail();
3518        ViewCell::NewMail();
3519        numViewCells = 0;
3520
3521        KdLeaf *leaf = tData.mNode;
3522
3523        RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
3524
3525        for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
3526        {
3527                VssRay *ray = (*rit).mRay;
3528
3529                int newViewCells;
3530
3531                // if hitpoint with one of the objects is inside this node, we
3532                // evaluate the volume of the view cells seen by this ray
3533                if (ray->mOriginObject && (GetLeaf(ray->mOrigin, ray->mOriginNode) == leaf))
3534                {
3535            vol += PrepareHeuristics(*ray, newViewCells);
3536                        numViewCells += newViewCells;
3537                }
3538
3539                // count double if both hit points are within the kd node
3540                if (ray->mTerminationObject && (GetLeaf(ray->mTermination, ray->mTerminationNode) == leaf))
3541                {
3542                        vol += PrepareHeuristics(*ray, newViewCells);
3543                        numViewCells += newViewCells;
3544                }
3545        }
3546
3547        //Debug << "vol: " << vol << endl;
3548        return vol;
3549}
3550
3551
3552void OspTree::EvalHeuristicsContribution(KdLeaf *leaf,
3553                                                                                 const SortableEntry &ci,
3554                                                                                 float &volLeft,
3555                                                                                 float &volRight,
3556                                                                                 int &pvsLeft,
3557                                                                                 int &pvsRight,
3558                                                                                 int &viewCellsLeft,
3559                                                                                 int &viewCellsRight)
3560{
3561        Intersectable *obj = ci.mObject;
3562        VssRay *ray = ci.mRay;
3563       
3564        switch (ci.mType)
3565        {
3566                // add reverse pvs to left side of split plane
3567                case SortableEntry::BOX_MIN:
3568                        ++ pvsLeft;
3569                        break;
3570                       
3571                case SortableEntry::BOX_MAX:
3572                        -- pvsRight;
3573                        break;
3574
3575                // compute volume contribution from view cells
3576                case SortableEntry::BOX_INTERSECT:
3577                                if ((ray->mOriginObject && (GetLeaf(ray->mOrigin, ray->mOriginNode) == leaf)) ||
3578                                        (ray->mTerminationObject && (GetLeaf(ray->mTermination, ray->mTerminationNode) == leaf)))
3579                                {
3580                                        EvalVolumeContribution(*ray, volLeft, volRight, viewCellsLeft, viewCellsRight);
3581                                }
3582                        break;
3583                default:
3584                        cout << "should not come here" << endl;
3585                        break;
3586        }
3587
3588        //cout << "vol left: " << volLeft << " vol right " << volRight << endl;
3589}
3590
3591
3592void OspTree::EvalVolumeContribution(const VssRay &ray,
3593                                                                         float &volLeft,
3594                                                                         float &volRight,
3595                                                                         int &viewCellsLeft,
3596                                                                         int &viewCellsRight)
3597{
3598        ViewCellContainer viewCells;
3599
3600        mVspTree->GetViewCells(ray, viewCells);
3601
3602        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
3603       
3604        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
3605        {
3606                // add view cells volume to left child node
3607                ViewCell *viewCell = *vit;
3608
3609                if (!viewCell->Mailed())
3610                {
3611                        viewCell->Mail();
3612                       
3613                        volLeft += viewCell->GetVolume();
3614
3615                        ++ viewCellsLeft;
3616                }
3617
3618                // remove from right child node
3619                if (-- viewCell->mCounter == 0)
3620                {
3621                        volRight -= viewCell->GetVolume();
3622
3623                        -- viewCellsRight;
3624                }
3625        }
3626}
3627
3628
3629void OspTree::SetViewCellsManager(ViewCellsManager *vcm)
3630{
3631        mViewCellsManager = vcm;
3632}
3633
3634
3635AxisAlignedBox3 OspTree::GetBoundingBox() const
3636{
3637        return mBoundingBox;
3638}
3639
3640
3641float OspTree::SelectSplitPlane(const OspTraversalData &tData,
3642                                                                AxisAlignedPlane &plane,
3643                                                                float &pFront,
3644                                                                float &pBack)
3645{
3646        float nPosition[3];
3647        float nCostRatio[3];
3648        float nProbFront[3];
3649        float nProbBack[3];
3650
3651        // create bounding box of node geometry
3652        AxisAlignedBox3 box = tData.mBoundingBox;
3653               
3654        int sAxis = 0;
3655        int bestAxis = -1;
3656
3657        if (mOnlyDrivingAxis)
3658        {
3659                sAxis = box.Size().DrivingAxis();
3660        }
3661
3662        // -- evaluate split cost for all three axis
3663        for (int axis = 0; axis < 3; ++ axis)
3664        {
3665                if (!mOnlyDrivingAxis || (axis == sAxis))
3666                {
3667                        if (1 || mUseCostHeuristics)
3668                        {
3669                                //-- place split plane using heuristics
3670                                int objectsFront, objectsBack;
3671
3672                                nCostRatio[axis] =
3673                                        EvalLocalCostHeuristics(tData,
3674                                                                           tData.mBoundingBox,
3675                                                                           axis,
3676                                                                           nPosition[axis],
3677                                                                           objectsFront,
3678                                                                           objectsBack);                       
3679                        }
3680                        /*      else
3681                        {
3682                                //-- split plane position is spatial median
3683
3684                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
3685
3686                                nCostRatio[axis] = EvalLocalSplitCost(tData,
3687                                                                                                          box,
3688                                                                                                          axis,
3689                                                                                                          nPosition[axis],
3690                                                                                                          nProbFront[axis],
3691                                                                                                          nProbBack[axis]);
3692                        }*/
3693                                               
3694                        if (bestAxis == -1)
3695                        {
3696                                bestAxis = axis;
3697                        }
3698                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
3699                        {
3700                                bestAxis = axis;
3701                        }
3702                }
3703        }
3704
3705        //-- assign values
3706
3707        plane.mAxis = bestAxis;
3708        // split plane position
3709    plane.mPosition = nPosition[bestAxis];
3710
3711        pFront = nProbFront[bestAxis];
3712        pBack = nProbBack[bestAxis];
3713
3714        Debug << "val: " << nCostRatio[bestAxis] << " axis: " << bestAxis << endl;
3715        return nCostRatio[bestAxis];
3716}
3717
3718
3719bool OspTree::EndPointInsideNode(KdLeaf *leaf,
3720                                                                 VssRay &ray,
3721                                                                 bool isTermination) const
3722{
3723        // test if the leaf where the hitpoint is located is the current leaf
3724        if (isTermination)
3725        {
3726                 return ray.mTerminationObject && (GetLeaf(ray.mTermination, ray.mOriginNode) == leaf);
3727        }
3728        else
3729        {
3730                 return ray.mOriginObject && (GetLeaf(ray.mOrigin, ray.mOriginNode) == leaf);
3731        }
3732}
3733
3734
3735#if TODO
3736void OspTree::EvalSubdivisionStats(const VspTraversalData &tData,
3737                                                                   const VspTraversalData &tFrontData,
3738                                                                   const VspTraversalData &tBackData
3739                                                                   )
3740{
3741        const float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
3742        const float cBack = (float)tBackData.mPvs * tBackData.mProbability;
3743        const float cData = (float)tData.mPvs * tData.mProbability;
3744       
3745        const float costDecr =
3746                (cFront + cBack - cData) / mBoundingBox.GetVolume();
3747
3748        mTotalCost += costDecr;
3749        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
3750
3751        AddSubdivisionStats(mVspStats.Leaves(),
3752                                                -costDecr,
3753                                                sc->GetPriority(),
3754                                                mTotalCost,
3755                                                (float)mTotalPvsSize / (float)mVspStats.Leaves());
3756}
3757#endif
3758
3759
3760int OspTree::ClassifyRay(VssRay *ray,
3761                                                 KdLeaf *leaf,
3762                                                 const AxisAlignedPlane &plane) const
3763{
3764        const bool originInside = EndPointInsideNode(leaf, *ray, false);
3765    const bool terminationInside = EndPointInsideNode(leaf, *ray, true);
3766
3767        const bool originGt =
3768                ray->mOrigin[plane.mAxis] > plane.mPosition;
3769               
3770        const bool terminationGt =
3771                ray->mTermination[plane.mAxis] > plane.mPosition;
3772
3773        // add view cell volume to front volume
3774        const bool addToFront = ((originInside && originGt) || (terminationInside && terminationGt));
3775
3776        // add view cell volume to back volume
3777        const bool addToBack = ((originInside && !originGt) || (terminationInside && !terminationGt));
3778
3779        // classify ray with respect to the child nodes the view cells contribute
3780        if (addToFront && addToBack)
3781                return 0;
3782        else if (addToBack)
3783                return -1;
3784        else if (addToFront)
3785                return 1;
3786
3787        return -2;
3788}
3789
3790
3791float OspTree::EvalRenderCostDecrease(const AxisAlignedPlane &candidatePlane,
3792                                                                          const OspTraversalData &tData,
3793                                                                          float &normalizedOldRenderCost) const
3794{
3795        float pvsFront = 0;
3796        float pvsBack = 0;
3797        float totalPvs = 0;
3798
3799        // probability that view point lies in back / front node
3800        float pOverall = 0;//data.mProbability; // todo matt§: value ok?
3801        float pFront = 0;
3802        float pBack = 0;
3803
3804        const float viewSpaceVol = mVspTree->GetBoundingBox().GetVolume();
3805
3806        //Intersectable::NewMail();
3807        KdLeaf::NewMail();
3808        ViewCell::NewMail();
3809
3810        KdLeaf *leaf = tData.mNode;
3811        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
3812
3813        // evaluate reverse pvs and view cell volume on left and right cell
3814        // note: should I take all leaf objects or rather the objects hit by rays?
3815        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
3816        {
3817                Intersectable *obj = *oit;
3818                const AxisAlignedBox3 box = obj->GetBox();
3819               
3820                ++ totalPvs;
3821
3822                //cout << "totalpvs " << totalPvs << endl;
3823
3824                // test if box falls in left / right child node
3825                if (box.Max(candidatePlane.mAxis) > candidatePlane.mPosition)
3826                {
3827                        ++ pvsFront;
3828                }
3829                if (box.Min(candidatePlane.mAxis) < candidatePlane.mPosition)
3830                {
3831                        ++ pvsBack;
3832                }
3833        }
3834
3835
3836        // sum up volume seen from the objects of left and right children
3837        // => the volume is the weight for the render cost equation
3838        ViewCell::NewMail(3);
3839       
3840        RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
3841
3842        for (rit = tData.mRays->begin(); rit < rit_end; ++ rit)
3843        {
3844                VssRay *ray = (*rit).mRay;
3845
3846                // test if intersection point with one of the objects is inside this node
3847                const bool originInside = EndPointInsideNode(leaf, *ray, false);
3848        const bool terminationInside = EndPointInsideNode(leaf, *ray, true);
3849
3850                if (originInside || terminationInside)
3851                {
3852                        // add volume to volumes of left and / or right children
3853                        // if one of the ray end points is inside
3854                        const int classification = ClassifyRay(ray, leaf, candidatePlane);
3855
3856                        ViewCellContainer viewCells;           
3857                        mVspTree->GetViewCells(*ray, viewCells);
3858                       
3859                        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
3860
3861                        // evaluate view cells volume contribution
3862                        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
3863                        {
3864                                AddViewCellVolume(*vit, classification, pFront, pBack, pOverall);
3865                        }
3866                }
3867        }
3868
3869
3870        //-- pvs rendering heuristics
3871        const float oldRenderCost = pOverall * totalPvs;
3872        const float newRenderCost = pvsFront * pFront + pvsBack * pBack;
3873
3874        // normalize volume with view space volume
3875        const float renderCostDecrease = (oldRenderCost - newRenderCost) / viewSpaceVol;
3876
3877        Debug << "\n==== eval render cost decrease ===" << endl
3878                  << "back pvs: " << pvsBack << " front pvs " << pvsFront << " total pvs: " << totalPvs << endl
3879                  << "back p: " << pBack / viewSpaceVol << " front p " << pFront / viewSpaceVol << " p: " << pOverall / viewSpaceVol << endl
3880                  << "old rc: " << oldRenderCost / viewSpaceVol << " new rc: " << newRenderCost / viewSpaceVol << endl
3881                  << "render cost decrease: " << renderCostDecrease << endl;
3882
3883        normalizedOldRenderCost = oldRenderCost / viewSpaceVol;
3884
3885        return renderCostDecrease;
3886}
3887
3888
3889void OspTree::ComputeBoundingBox(const ObjectContainer &objects,
3890                                                                 AxisAlignedBox3 *forcedBoundingBox)
3891{
3892        if (forcedBoundingBox)
3893        {
3894                mBoundingBox = *forcedBoundingBox;
3895        }
3896        else // compute vsp tree bounding box
3897        {
3898                mBoundingBox.Initialize();
3899
3900                ObjectContainer::const_iterator oit, oit_end = objects.end();
3901
3902                //-- compute bounding box
3903        for (oit = objects.begin(); oit != oit_end; ++ oit)
3904                {
3905                        Intersectable *obj = *oit;
3906
3907                        // compute bounding box of view space
3908                        mBoundingBox.Include(obj->GetBox());
3909                        mBoundingBox.Include(obj->GetBox());
3910                }
3911        }
3912}
3913
3914
3915void OspTree::ProcessMultipleRefs(KdLeaf *leaf) const
3916{
3917        // find objects from multiple kd-leaves
3918        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
3919
3920        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
3921        {
3922                Intersectable *object = *oit;
3923               
3924                if (object->mReferences > 1)
3925                {
3926                        leaf->mMultipleObjects.push_back(object);
3927                }
3928        }
3929}
3930
3931
3932void OspTree::CollectLeaves(vector<KdLeaf *> &leaves) const
3933{
3934        stack<KdNode *> nodeStack;
3935        nodeStack.push(mRoot);
3936
3937        while (!nodeStack.empty())
3938        {
3939                KdNode *node = nodeStack.top();
3940                nodeStack.pop();
3941                if (node->IsLeaf())
3942                {
3943                        KdLeaf *leaf = (KdLeaf *)node;
3944                        leaves.push_back(leaf);
3945                }
3946                else
3947                {
3948                        KdInterior *interior = (KdInterior *)node;
3949                        nodeStack.push(interior->mBack);
3950                        nodeStack.push(interior->mFront);
3951                }
3952        }
3953}
3954
3955
3956AxisAlignedBox3 OspTree::GetBoundingBox(KdNode *node) const
3957{
3958        if (!node->mParent)
3959                return mBoundingBox;
3960
3961        if (!node->IsLeaf())
3962        {
3963                return (dynamic_cast<KdInterior *>(node))->mBox;               
3964        }
3965
3966        KdInterior *parent = dynamic_cast<KdInterior *>(node->mParent);
3967
3968        AxisAlignedBox3 box(parent->mBox);
3969
3970        if (parent->mFront == node)
3971                box.SetMin(parent->mAxis, parent->mPosition);
3972    else
3973                box.SetMax(parent->mAxis, parent->mPosition);
3974
3975        return box;
3976}
3977
3978
3979void OspTree::CollectViewCells(KdLeaf *leaf, ViewCellContainer &viewCells)
3980{
3981        ObjectContainer::const_iterator oit, oit_end = leaf->mObjects.end();
3982
3983        ViewCell::NewMail();
3984
3985        // loop through all object and collect view cell pvs of this node
3986        for (oit = leaf->mObjects.begin(); oit != oit_end; ++ oit)
3987        {
3988                Intersectable *obj = *oit;
3989
3990                ViewCellPvsMap::const_iterator vit, vit_end = obj->mViewCellPvs.mEntries.end();
3991
3992                for (vit = obj->mViewCellPvs.mEntries.begin(); vit != vit_end; ++ vit)
3993                {
3994            ViewCell *vc = (*vit).first;
3995
3996                        if (!vc->Mailed())
3997                        {
3998                                vc->Mail();
3999                                viewCells.push_back(vc);
4000                        }
4001                }
4002        }
4003}
4004
4005
4006void OspTree::CollectDirtyCandidates(OspSplitCandidate *sc,
4007                                                                         vector<SplitCandidate *> &dirtyList)
4008{
4009        OspTraversalData &tData = sc->mParentData;
4010        KdLeaf *node = tData.mNode;
4011       
4012        ViewCell::NewMail();
4013        ViewCellContainer viewCells;
4014        ViewCellContainer tmpViewCells;
4015
4016        RayInfoContainer::const_iterator rit, rit_end = tData.mRays->end();
4017
4018        // find all view cells associated with the rays, add them to view cells
4019        for (rit = tData.mRays->begin(); rit != rit_end; ++ rit)
4020        {
4021                VssRay *ray = (*rit).mRay;
4022                mVspTree->GetViewCells(*ray, tmpViewCells);
4023
4024                ViewCellContainer::const_iterator vit, vit_end = tmpViewCells.end();
4025
4026                for (vit = tmpViewCells.begin(); vit != vit_end; ++ vit)
4027                {
4028                        VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);
4029
4030                        if (!vc->Mailed())
4031                        {
4032                                vc->Mail();
4033                                viewCells.push_back(vc);
4034                        }
4035                }
4036        }
4037
4038
4039        // split candidates holding this view cells
4040        // are thrown into dirty list
4041        ViewCellContainer::const_iterator vit, vit_end = viewCells.end();
4042
4043        for (vit = viewCells.begin(); vit != vit_end; ++ vit)
4044        {
4045                VspViewCell *vc = dynamic_cast<VspViewCell *>(*vit);
4046
4047                VspLeaf *leaf = vc->mLeaf;
4048                dirtyList.push_back(leaf->GetSplitCandidate());
4049        }
4050}
4051
4052
4053int OspTree::ComputePvsSize(const ObjectContainer &objects) const
4054{
4055        int pvsSize = 0;
4056
4057        ObjectContainer::const_iterator oit, oit_end = objects.end();
4058
4059        for (oit = objects.begin(); oit != oit_end; ++ oit)
4060        {
4061                pvsSize += (*oit)->mViewCellPvs.GetSize();
4062        }
4063
4064        return pvsSize;
4065}
4066
4067
4068KdNode *OspTree::GetRoot() const
4069{
4070        return mRoot;
4071}
4072
4073
4074KdLeaf *OspTree::GetLeaf(const Vector3 &pt, KdNode *node) const
4075{
4076        // start from root of tree
4077        if (node == NULL)
4078        {
4079                node = mRoot;
4080        }
4081
4082        stack<KdNode *> nodeStack;
4083        nodeStack.push(node);
4084 
4085        KdLeaf *leaf = NULL;
4086 
4087        while (!nodeStack.empty()) 
4088        {
4089                KdNode *node = nodeStack.top();
4090                nodeStack.pop();
4091       
4092                if (node->IsLeaf())
4093                {
4094                        leaf = dynamic_cast<KdLeaf *>(node);
4095                }
4096                else   
4097                {       
4098                        // find point
4099                        KdInterior *interior = dynamic_cast<KdInterior *>(node);
4100       
4101                        if (interior->mPosition > pt[interior->mAxis])
4102                        {
4103                                nodeStack.push(interior->mBack);
4104                        }
4105                        else
4106                        {
4107                                nodeStack.push(interior->mFront);
4108                        }
4109                }
4110        }
4111 
4112        return leaf;
4113}
4114
4115
4116void OspTree::ProcessRays(const VssRayContainer &sampleRays,
4117                                                  RayInfoContainer &rays)
4118{
4119        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
4120
4121        long startTime = GetTime();
4122
4123        cout << "storing rays ... ";
4124
4125        Intersectable::NewMail();
4126
4127        //-- store rays and objects
4128        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
4129        {
4130                VssRay *ray = *rit;
4131                float minT, maxT;
4132                static Ray hray;
4133
4134                hray.Init(*ray);
4135               
4136                // TODO: not very efficient to implictly cast between rays types
4137                if (GetBoundingBox().GetRaySegment(hray, minT, maxT))
4138                {
4139                        float len = ray->Length();
4140
4141                        if (!len)
4142                                len = Limits::Small;
4143
4144                        rays.push_back(RayInfo(ray, minT / len, maxT / len));
4145                }
4146        }
4147
4148        cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
4149}
4150
4151
4152
4153int OspTree::SplitRays(const AxisAlignedPlane &plane,
4154                                           RayInfoContainer &rays,
4155                                           RayInfoContainer &frontRays,
4156                                           RayInfoContainer &backRays) const
4157{
4158        int splits = 0;
4159
4160        RayInfoContainer::const_iterator rit, rit_end = rays.end();
4161
4162        for (rit = rays.begin(); rit != rit_end; ++ rit)
4163        {
4164                RayInfo bRay = *rit;
4165               
4166                VssRay *ray = bRay.mRay;
4167                float t;
4168
4169                // get classification and receive new t
4170                //-- test if start point behind or in front of plane
4171                const int side = bRay.ComputeRayIntersection(plane.mAxis, plane.mPosition, t);
4172                       
4173                if (side == 0)
4174                {
4175                        ++ splits;
4176
4177                        if (ray->HasPosDir(plane.mAxis))
4178                        {
4179                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
4180                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
4181                        }
4182                        else
4183                        {
4184                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
4185                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
4186                        }
4187                }
4188                else if (side == 1)
4189                {
4190                        frontRays.push_back(bRay);
4191                }
4192                else
4193                {
4194                        backRays.push_back(bRay);
4195                }
4196        }
4197
4198        return splits;
4199}
4200
4201
4202KdIntersectable *OspTree::GetOrCreateKdIntersectable(KdNode *node)
4203{
4204        // search nodes
4205        std::map<KdNode *, KdIntersectable *>::
4206                const_iterator it = mKdIntersectables.find(node);
4207
4208        if (it != mKdIntersectables.end())
4209        {
4210                return (*it).second;
4211        }
4212
4213        // not in map => create new entry
4214        KdIntersectable *kdObj= new KdIntersectable(node);
4215
4216        mKdIntersectables[node] = kdObj;
4217
4218        return kdObj;
4219}
4220
4221
4222void OspTree::AddViewCellVolume(ViewCell *vc,
4223                                                                const int cf,
4224                                                                float &frontVol,
4225                                                                float &backVol,
4226                                                                float &totalVol) const
4227{
4228        //const float renderCost = mViewCellsManager->SimpleRay &raynderCost(obj);
4229        const float vol = vc->GetVolume();
4230
4231        // view cell not found yet => new
4232        if (!vc->Mailed() && !vc->Mailed(1) && !vc->Mailed(2))
4233        {
4234                totalVol += vol;
4235        }
4236
4237        if (cf >= 0) // front volume
4238        {
4239                if (!vc->Mailed() && !vc->Mailed(2))
4240                {
4241                        frontVol += vol;
4242               
4243                        // already in back volume => in both volumes
4244                        if (vc->Mailed(1))
4245                                vc->Mail(2);
4246                        else
4247                                vc->Mail();
4248                }
4249        }
4250
4251        if (cf <= 0) // back volume
4252        {
4253                if (!vc->Mailed(1) && !vc->Mailed(2))
4254                {
4255                        backVol += vol;
4256               
4257                        // already in front volume => in both volume
4258                        if (vc->Mailed())
4259                                vc->Mail(2);
4260                        else
4261                                vc->Mail(1);
4262                }
4263        }
4264}
4265
4266
4267/********************************************************************/
4268/*               class HierarchyManager implementation              */
4269/********************************************************************/
4270
4271
4272HierarchyManager::HierarchyManager(VspTree &vspTree, OspTree &ospTree):
4273mVspTree(vspTree), mOspTree(ospTree)
4274{
4275        // cross references
4276        mVspTree.mOspTree = &ospTree;
4277        mOspTree.mVspTree = &vspTree;
4278}
4279
4280
4281SplitCandidate *HierarchyManager::NextSplitCandidate()
4282{
4283        SplitCandidate *splitCandidate = mTQueue.Top();
4284
4285        cout << "next candidate: " << splitCandidate->Type()
4286                << ", priority: " << splitCandidate->GetPriority() << endl;
4287
4288        mTQueue.Pop();
4289
4290        return splitCandidate;
4291}
4292
4293
4294VspTree::VspSplitCandidate *HierarchyManager::PrepareVsp(const VssRayContainer &sampleRays,
4295                                                                                         AxisAlignedBox3 *forcedViewSpace,
4296                                                                                         RayInfoContainer &rays)
4297{
4298        // store pointer to this tree
4299        VspTree::VspSplitCandidate::sVspTree = &mVspTree;
4300        mVspTree.mVspStats.nodes = 1;
4301
4302        // compute view space bounding box
4303        mVspTree.ComputeBoundingBox(sampleRays, forcedViewSpace);
4304
4305        // initialise termination criteria
4306        mVspTree.mTermMinProbability *= mBoundingBox.GetVolume();
4307        mVspTree.mGlobalCostMisses = 0;
4308
4309        // get clipped rays
4310        mVspTree.ProcessRays(sampleRays, rays);
4311
4312        const int pvsSize = mVspTree.ComputePvsSize(rays);
4313       
4314        Debug <<  "pvs size: " << (int)pvsSize << endl;
4315        Debug <<  "rays size: " << (int)rays.size() << endl;
4316
4317        //-- prepare view space partition
4318
4319        // add first candidate for view space partition
4320        VspLeaf *leaf = new VspLeaf();
4321        mVspTree.mRoot = leaf;
4322
4323        const float prop = mVspTree.mBoundingBox.GetVolume();
4324
4325        // first vsp traversal data
4326        VspTree::VspTraversalData vData(leaf,
4327                                                                        0,
4328                                                                        &rays,
4329                                                                        pvsSize,       
4330                                                                        prop,
4331                                                                        mVspTree.mBoundingBox);
4332
4333
4334        // create first view cell
4335        mVspTree.CreateViewCell(vData, false);
4336               
4337        // add first view cell to all the objects view cell pvs
4338        ObjectPvsMap::const_iterator oit,
4339                oit_end = leaf->GetViewCell()->GetPvs().mEntries.end();
4340
4341
4342        for (oit = leaf->GetViewCell()->GetPvs().mEntries.begin(); oit != oit_end; ++ oit)
4343        {
4344                Intersectable *obj = (*oit).first;
4345                obj->mViewCellPvs.AddSample(leaf->GetViewCell(), 1);
4346        }
4347
4348        // compute first split candidate
4349        VspTree::VspSplitCandidate *splitCandidate =
4350                new VspTree::VspSplitCandidate(vData);
4351    mVspTree.EvalSplitCandidate(*splitCandidate);
4352
4353        return splitCandidate;
4354}
4355
4356
4357OspTree::OspSplitCandidate * HierarchyManager::PrepareOsp(const VssRayContainer &sampleRays,
4358                                                                                          const ObjectContainer &objects,
4359                                                                                          AxisAlignedBox3 *forcedObjectSpace,
4360                                                                                          RayInfoContainer &rays)
4361{
4362        // store pointer to this tree
4363        OspTree::OspSplitCandidate::sOspTree = &mOspTree;
4364        mOspTree.mOspStats.nodes = 1;
4365       
4366        // compute bounding box from objects
4367        mOspTree.ComputeBoundingBox(objects, forcedObjectSpace);
4368
4369        mOspTree.mTermMinProbability *= mBoundingBox.GetVolume();
4370        mGlobalCostMisses = 0;
4371
4372        // get clipped rays
4373        mOspTree.ProcessRays(sampleRays, rays);
4374
4375        // add first candidate for view space partition
4376        KdLeaf *kdleaf = new KdLeaf(NULL, 0);
4377        kdleaf->mObjects = objects;
4378
4379        mOspTree.mRoot = kdleaf;
4380       
4381        // TODO matt: this is something different than pvs size.
4382        const int pvsSize = mOspTree.ComputePvsSize(objects);
4383        const float prop = mOspTree.mBoundingBox.GetVolume();
4384
4385        // first osp traversal data
4386        OspTree::OspTraversalData oData(kdleaf,
4387                                                                        0,
4388                                                                        &rays,
4389                                                                        pvsSize,
4390                                                                        prop,
4391                                                                        mOspTree.mBoundingBox);
4392
4393               
4394        // compute first split candidate
4395        OspTree::OspSplitCandidate *oSplitCandidate =
4396                new OspTree::OspSplitCandidate(oData);
4397
4398    mOspTree.EvalSplitCandidate(*oSplitCandidate);
4399
4400        return oSplitCandidate;
4401}
4402
4403
4404void HierarchyManager::PrepareConstruction(const VssRayContainer &sampleRays,
4405                                                                                   const ObjectContainer &objects,
4406                                                                                   AxisAlignedBox3 *forcedViewSpace,
4407                                                                                   RayInfoContainer &viewSpaceRays,
4408                                                                                   RayInfoContainer &objectSpaceRays)
4409{
4410        SplitCandidate *vsc =
4411                PrepareVsp(sampleRays, forcedViewSpace, viewSpaceRays);
4412        mTQueue.Push(vsc);
4413
4414        SplitCandidate *osc =
4415                PrepareOsp(sampleRays, objects, forcedViewSpace, objectSpaceRays);
4416
4417        mTQueue.Push(osc);
4418}
4419
4420
4421bool HierarchyManager::GlobalTerminationCriteriaMet(SplitCandidate *candidate) const
4422{
4423        // TODO matt: make virtual by creating superclasss for traveral data
4424        return candidate->GlobalTerminationCriteriaMet();
4425}
4426
4427
4428void HierarchyManager::Construct(const VssRayContainer &sampleRays,
4429                                                                 const ObjectContainer &objects,
4430                                                                 AxisAlignedBox3 *forcedViewSpace)
4431{
4432        RayInfoContainer *objectSpaceRays = new RayInfoContainer();
4433        RayInfoContainer *viewSpaceRays = new RayInfoContainer();
4434
4435        // prepare vsp and osp trees for traversal
4436        PrepareConstruction(sampleRays, objects, forcedViewSpace, *viewSpaceRays, *objectSpaceRays);
4437
4438        mVspTree.mVspStats.Reset();
4439        mVspTree.mVspStats.Start();
4440
4441        cout << "Constructing view space / object space tree ... \n";
4442        const long startTime = GetTime();       
4443       
4444        RunConstruction(true);
4445
4446        cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl;
4447
4448        mVspTree.mVspStats.Stop();
4449}
4450
4451
4452bool HierarchyManager::SubdivideSplitCandidate(SplitCandidate *sc)
4453{
4454        const bool globalTerminationCriteriaMet =
4455                        GlobalTerminationCriteriaMet(sc);
4456
4457        const bool vspSplit = (sc->Type() == SplitCandidate::VIEW_SPACE);
4458
4459        if (vspSplit)
4460        {
4461                mVspTree.Subdivide(mTQueue, sc, globalTerminationCriteriaMet);
4462                return true;
4463        }
4464        else
4465        {
4466                mOspTree.Subdivide(mTQueue, sc, globalTerminationCriteriaMet);
4467                return false;
4468        }
4469}
4470
4471
4472void HierarchyManager::RunConstruction(const bool repair)
4473{
4474        int numNodes = 0;
4475
4476        while (!FinishedConstruction())
4477        {
4478                SplitCandidate *splitCandidate = NextSplitCandidate();
4479           
4480                cout << "nodes: " << ++ numNodes << endl;
4481
4482                mTotalCost -= splitCandidate->GetRenderCostDecrease();
4483
4484                //-- subdivide leaf node
4485                SubdivideSplitCandidate(splitCandidate);
4486               
4487                // cost ratio of cost decrease / totalCost
4488                const float costRatio = splitCandidate->GetRenderCostDecrease() / mTotalCost;
4489               
4490                Debug << "\n**********" << endl
4491                          << "total cost: " << mTotalCost << " render cost decr: "
4492                          << splitCandidate->GetRenderCostDecrease()
4493                          << " cost ratio: " << costRatio << endl << endl;
4494
4495                if (costRatio < mTermMinGlobalCostRatio)
4496                        ++ mGlobalCostMisses;
4497
4498                // reevaluate candidates affected by the split
4499                // for view space splits, this would be object space splits
4500                // and other way round
4501                if (repair)
4502                        RepairQueue();
4503
4504                DEL_PTR(splitCandidate);
4505        }
4506}
4507
4508
4509void HierarchyManager::Construct2(const VssRayContainer &sampleRays,
4510                                                                  const ObjectContainer &objects,
4511                                                                  AxisAlignedBox3 *forcedViewSpace)
4512{
4513        RayInfoContainer *viewSpaceRays = new RayInfoContainer();
4514        RayInfoContainer *objectSpaceRays = new RayInfoContainer();
4515
4516
4517        /////////////////////////////////////////////////////////////
4518        // view space space partition
4519        /////////////////////////////////////////////////////////////
4520
4521
4522        // makes no sense otherwise because only one kd cell available
4523        // during view space partition
4524        const bool savedCountMethod = mVspTree.mUseKdPvsForHeuristics;
4525        const bool savedStoreMethod = mVspTree.mStoreKdPvs;
4526       
4527        mVspTree.mUseKdPvsForHeuristics = false;
4528        mVspTree.mStoreKdPvs = false;
4529
4530        VspTree::VspSplitCandidate *vsc =
4531                PrepareVsp(sampleRays, forcedViewSpace, *viewSpaceRays);
4532        // add to queue
4533        mTQueue.Push(vsc);
4534
4535        long startTime = GetTime();
4536        cout << "starting vsp contruction ... " << endl;
4537
4538        mVspTree.mVspStats.Reset();
4539        mVspTree.mVspStats.Start();
4540
4541        int i = 0;
4542
4543        // all objects can be seen from everywhere
4544        mTotalCost = (float)vsc->mParentData.mPvs;
4545
4546        const bool repairQueue = false;
4547
4548        // process view space candidates
4549        RunConstruction(repairQueue);
4550
4551        cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl;
4552        mVspTree.mVspStats.Stop();
4553       
4554
4555
4556        /////////////////////////////////////////////////////////////
4557        // object space partition
4558        /////////////////////////////////////////////////////////////
4559
4560        Debug << "\n$$$$$$$$$ osp tree construction $$$$$$$$$$\n" << endl;
4561        cout << "starting osp contruction ... " << endl;
4562
4563        // compute first candidate
4564        OspTree::OspSplitCandidate *osc =
4565                PrepareOsp(sampleRays, objects, forcedViewSpace, *objectSpaceRays);
4566
4567    mTQueue.Push(osc);
4568
4569        mOspTree.mOspStats.Reset();
4570        mOspTree.mOspStats.Start();
4571
4572        startTime = GetTime();
4573
4574        // reset cost
4575        // start with one big kd cell - all objects can be seen from everywhere
4576        // note: only true for view space = object space
4577        mTotalCost = (float)osc->mParentData.mNode->mObjects.size();//(float)sc->mParentData.mPvs;
4578
4579        Debug << "reseting cost, new total cost: " << mTotalCost << endl;
4580
4581        // process object space candidates
4582        RunConstruction(repairQueue);
4583       
4584        cout << "finished in " << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
4585
4586        mOspTree.mOspStats.Stop();
4587
4588        // reset parameters
4589        mVspTree.mUseKdPvsForHeuristics = savedCountMethod;
4590        mVspTree.mStoreKdPvs = savedStoreMethod;
4591}
4592
4593
4594void HierarchyManager::Construct3(const VssRayContainer &sampleRays,
4595                                                                  const ObjectContainer &objects,
4596                                                                  AxisAlignedBox3 *forcedViewSpace)
4597{
4598        // only view space partition
4599        // object kd tree is taken for osp
4600
4601        mVspTree.mVspStats.Reset();
4602        mVspTree.mVspStats.Start();
4603
4604        RayInfoContainer *viewSpaceRays = new RayInfoContainer();
4605       
4606        SplitCandidate *sc = PrepareVsp(sampleRays, forcedViewSpace, *viewSpaceRays);
4607        mTQueue.Push(sc);
4608
4609        cout << "starting vsp contruction ... " << endl;
4610
4611        long startTime = GetTime();
4612
4613        RunConstruction(false);
4614
4615        cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl;
4616        mVspTree.mVspStats.Stop();
4617}
4618
4619
4620bool HierarchyManager::FinishedConstruction() const
4621{
4622        return mTQueue.Empty();
4623}
4624
4625
4626void HierarchyManager::CollectDirtyCandidates(vector<SplitCandidate *> &dirtyList)
4627{       
4628        // we have either a object space or view space split
4629        if (mCurrentCandidate->Type() == SplitCandidate::VIEW_SPACE)
4630        {
4631                VspTree::VspSplitCandidate *sc =
4632                        dynamic_cast<VspTree::VspSplitCandidate *>(mCurrentCandidate);
4633
4634                mVspTree.CollectDirtyCandidates(sc, dirtyList);
4635        }
4636        else // object space split
4637        {                       
4638                OspTree::OspSplitCandidate *sc =
4639                        dynamic_cast<OspTree::OspSplitCandidate *>(mCurrentCandidate);
4640
4641                mOspTree.CollectDirtyCandidates(sc, dirtyList);
4642        }
4643}
4644
4645
4646void HierarchyManager::RepairQueue()
4647{
4648        // TODO
4649        // for each update of the view space partition:
4650        // the candidates from object space partition which
4651        // have been afected by the view space split (the kd split candidates
4652        // which saw the view cell which was split) must be reevaluated
4653        // (maybe not locally, just reinsert them into the queue)
4654        //
4655        // vice versa for the view cells
4656        // for each update of the object space partition
4657        // reevaluate split candidate for view cells which saw the split kd cell
4658        //
4659        // the priority queue update can be solved by implementing a binary heap
4660        // (explicit data structure, binary tree)
4661        // *) inserting and removal is efficient
4662        // *) search is not efficient => store queue position with each
4663        // split candidate
4664
4665        // collect list of "dirty" candidates
4666        vector<SplitCandidate *> dirtyList;
4667        CollectDirtyCandidates(dirtyList);
4668       
4669        //-- reevaluate the dirty list
4670        vector<SplitCandidate *>::const_iterator sit, sit_end = dirtyList.end();
4671
4672        for (sit = dirtyList.begin(); sit != sit_end; ++ sit)
4673        {
4674                SplitCandidate* sc = *sit;
4675
4676                mTQueue.Erase(sc);
4677
4678                // reevaluate
4679                sc->EvalPriority();
4680
4681                // reinsert
4682                mTQueue.Push(sc);
4683        }
4684}
4685
4686}
Note: See TracBrowser for help on using the repository browser.