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

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