source: GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.cpp @ 729

Revision 729, 90.3 KB checked in by mattausch, 18 years ago (diff)

added cost flexible render cost measurements

Line 
1#include <stack>
2#include <time.h>
3#include <iomanip>
4
5#include "Plane3.h"
6#include "VspBspTree.h"
7#include "Mesh.h"
8#include "common.h"
9#include "ViewCell.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 "ViewCellBsp.h"
17#include "ViewCellsManager.h"
18#include "Beam.h"
19
20#define USE_FIXEDPOINT_T 0
21
22
23//-- static members
24
25
26int VspBspTree::sFrontId = 0;
27int VspBspTree::sBackId = 0;
28int VspBspTree::sFrontAndBackId = 0;
29
30typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;
31
32
33// pvs penalty can be different from pvs size
34inline float EvalPvsPenalty(const int pvs,
35                                                        const int lower,
36                                                        const int upper)
37{
38        // clamp to minmax values
39        if (pvs < lower)
40                return (float)lower;
41        if (pvs > upper)
42                return (float)upper;
43
44        return (float)pvs;
45}
46
47
48
49
50/******************************************************************************/
51/*                       class VspBspTree implementation                      */
52/******************************************************************************/
53
54
55VspBspTree::VspBspTree():
56mRoot(NULL),
57mUseAreaForPvs(false),
58mCostNormalizer(Limits::Small),
59mViewCellsManager(NULL),
60mOutOfBoundsCell(NULL),
61mStoreRays(false),
62mRenderCostWeight(0.5),
63mUseRandomAxis(false),
64mTimeStamp(1)
65{
66        bool randomize = false;
67        environment->GetBoolValue("VspBspTree.Construction.randomize", randomize);
68        if (randomize)
69                Randomize(); // initialise random generator for heuristics
70
71        //-- termination criteria for autopartition
72        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
73        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
74        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
75        environment->GetFloatValue("VspBspTree.Termination.minProbability", mTermMinProbability);
76        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
77        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
78        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
79        environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
80        environment->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
81
82        //-- max cost ratio for early tree termination
83        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
84
85        environment->GetFloatValue("VspBspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
86        environment->GetIntValue("VspBspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
87
88        // HACK//mTermMinPolygons = 25;
89
90        //-- factors for bsp tree split plane heuristics
91        environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
92        environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
93
94
95        //-- partition criteria
96        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
97        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
98        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
99
100        environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
101        environment->GetIntValue("VspBspTree.maxTests", mMaxTests);
102
103        // if only the driving axis is used for axis aligned split
104        environment->GetBoolValue("VspBspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
105       
106        //-- termination criteria for axis aligned split
107        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxRayContribution",
108                                                                mTermMaxRayContriForAxisAligned);
109        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
110                                                         mTermMinRaysForAxisAligned);
111
112        //environment->GetFloatValue("VspBspTree.maxTotalMemory", mMaxTotalMemory);
113        environment->GetFloatValue("VspBspTree.maxStaticMemory", mMaxMemory);
114
115        environment->GetFloatValue("VspBspTree.Construction.renderCostWeight", mRenderCostWeight);
116        environment->GetBoolValue("VspBspTree.usePolygonSplitIfAvailable", mUsePolygonSplitIfAvailable);
117
118        environment->GetBoolValue("VspBspTree.useCostHeuristics", mUseCostHeuristics);
119        environment->GetBoolValue("VspBspTree.useSplitCostQueue", mUseSplitCostQueue);
120        environment->GetBoolValue("VspBspTree.simulateOctree", mSimulateOctree);
121        environment->GetBoolValue("VspBspTree.useRandomAxis", mUseRandomAxis);
122        environment->GetBoolValue("VspBspTree.useBreathFirstSplits", mBreathFirstSplits);
123
124        environment->GetBoolValue("ViewCells.PostProcess.emptyViewCellsMerge", mEmptyViewCellsMergeAllowed);
125       
126        char subdivisionStatsLog[100];
127        environment->GetStringValue("VspBspTree.subdivisionStats", subdivisionStatsLog);
128        mSubdivisionStats.open(subdivisionStatsLog);
129
130        //-- debug output
131
132        Debug << "******* VSP BSP options ******** " << endl;
133    Debug << "max depth: " << mTermMaxDepth << endl;
134        Debug << "min PVS: " << mTermMinPvs << endl;
135        Debug << "min probabiliy: " << mTermMinProbability << endl;
136        Debug << "min rays: " << mTermMinRays << endl;
137        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
138        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
139        Debug << "miss tolerance: " << mTermMissTolerance << endl;
140        Debug << "max view cells: " << mMaxViewCells << endl;
141        Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
142        //Debug << "max plane candidates: " << mMaxRayCandidates << endl;
143        Debug << "randomize: " << randomize << endl;
144
145        Debug << "using area for pvs: " << mUseAreaForPvs << endl;
146        Debug << "render cost weight: " << mRenderCostWeight << endl;
147        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
148        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
149        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
150        Debug << "max memory: " << mMaxMemory << endl;
151        Debug << "use poly split if available: " << mUsePolygonSplitIfAvailable << endl;
152        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
153        Debug << "use split cost queue: " << mUseSplitCostQueue << endl;
154        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
155        Debug << "use random axis: " << mUseRandomAxis << endl;
156        Debug << "breath first splits: " << mBreathFirstSplits << endl;
157        Debug << "empty view cells merge: " << mEmptyViewCellsMergeAllowed << endl;
158
159        Debug << "octree: " << mSimulateOctree << endl;
160
161        Debug << "Split plane strategy: ";
162
163        if (mSplitPlaneStrategy & RANDOM_POLYGON)
164        {
165                Debug << "random polygon ";
166        }
167        if (mSplitPlaneStrategy & AXIS_ALIGNED)
168        {
169                Debug << "axis aligned ";
170        }
171        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
172        {
173                mCostNormalizer += mLeastRaySplitsFactor;
174                Debug << "least ray splits ";
175        }
176        if (mSplitPlaneStrategy & BALANCED_RAYS)
177        {
178                mCostNormalizer += mBalancedRaysFactor;
179                Debug << "balanced rays ";
180        }
181        if (mSplitPlaneStrategy & PVS)
182        {
183                mCostNormalizer += mPvsFactor;
184                Debug << "pvs";
185        }
186
187
188        mSplitCandidates = new vector<SortableEntry>;
189
190        Debug << endl;
191}
192
193
194BspViewCell *VspBspTree::GetOutOfBoundsCell()
195{
196        return mOutOfBoundsCell;
197}
198
199
200BspViewCell *VspBspTree::GetOrCreateOutOfBoundsCell()
201{
202        if (!mOutOfBoundsCell)
203        {
204                mOutOfBoundsCell = new BspViewCell();
205                mOutOfBoundsCell->SetId(-1);
206                mOutOfBoundsCell->SetValid(false);
207        }
208
209        return mOutOfBoundsCell;
210}
211
212
213const BspTreeStatistics &VspBspTree::GetStatistics() const
214{
215        return mBspStats;
216}
217
218
219VspBspTree::~VspBspTree()
220{
221        DEL_PTR(mRoot);
222        DEL_PTR(mSplitCandidates);
223}
224
225
226int VspBspTree::AddMeshToPolygons(Mesh *mesh,
227                                                                  PolygonContainer &polys,
228                                                                  MeshInstance *parent)
229{
230        FaceContainer::const_iterator fi;
231
232        // copy the face data to polygons
233        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
234        {
235                Polygon3 *poly = new Polygon3((*fi), mesh);
236
237                if (poly->Valid(mEpsilon))
238                {
239                        poly->mParent = parent; // set parent intersectable
240                        polys.push_back(poly);
241                }
242                else
243                        DEL_PTR(poly);
244        }
245        return (int)mesh->mFaces.size();
246}
247
248
249int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells,
250                                                          PolygonContainer &polys,
251                                                          int maxObjects)
252{
253        int limit = (maxObjects > 0) ?
254                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
255
256        int polysSize = 0;
257
258        for (int i = 0; i < limit; ++ i)
259        {
260                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
261                {
262                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
263                        polysSize +=
264                                AddMeshToPolygons(viewCells[i]->GetMesh(),
265                                                                  polys,
266                                                                  viewCells[i]);
267                }
268        }
269
270        return polysSize;
271}
272
273
274int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects,
275                                                                 PolygonContainer &polys,
276                                                                 int maxObjects)
277{
278        int limit = (maxObjects > 0) ?
279                Min((int)objects.size(), maxObjects) : (int)objects.size();
280
281        for (int i = 0; i < limit; ++i)
282        {
283                Intersectable *object = objects[i];//*it;
284                Mesh *mesh = NULL;
285
286                switch (object->Type()) // extract the meshes
287                {
288                case Intersectable::MESH_INSTANCE:
289                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
290                        break;
291                case Intersectable::VIEW_CELL:
292                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
293                        break;
294                        // TODO: handle transformed mesh instances
295                default:
296                        Debug << "intersectable type not supported" << endl;
297                        break;
298                }
299
300        if (mesh) // copy the mesh data to polygons
301                {
302                        mBox.Include(object->GetBox()); // add to BSP tree aabb
303                        AddMeshToPolygons(mesh, polys, NULL);
304                }
305        }
306
307        return (int)polys.size();
308}
309
310
311void VspBspTree::Construct(const VssRayContainer &sampleRays,
312                                                   AxisAlignedBox3 *forcedBoundingBox)
313{
314        mBspStats.nodes = 1;
315        mBox.Initialize();      // initialise BSP tree bounding box
316
317        if (forcedBoundingBox)
318                mBox = *forcedBoundingBox;
319       
320        PolygonContainer polys;
321        RayInfoContainer *rays = new RayInfoContainer();
322
323        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
324
325        long startTime = GetTime();
326
327        cout << "Extracting polygons from rays ... ";
328
329        Intersectable::NewMail();
330
331        int numObj = 0;
332
333        //-- extract polygons intersected by the rays
334        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
335        {
336                VssRay *ray = *rit;
337
338                if ((mBox.IsInside(ray->mTermination) || !forcedBoundingBox) &&
339                        ray->mTerminationObject &&
340                        !ray->mTerminationObject->Mailed())
341                {
342                        ray->mTerminationObject->Mail();
343                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
344                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
345                        ++ numObj;
346
347                        //-- compute bounding box
348                        if (!forcedBoundingBox)
349                                mBox.Include(ray->mTermination);
350                }
351
352                if ((mBox.IsInside(ray->mOrigin) || !forcedBoundingBox) &&
353                        ray->mOriginObject &&
354                        !ray->mOriginObject->Mailed())
355                {
356                        ray->mOriginObject->Mail();
357                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
358                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
359                        ++ numObj;
360
361                        //-- compute bounding box
362                        if (!forcedBoundingBox)
363                                mBox.Include(ray->mOrigin);
364                }
365        }
366       
367        Debug << "maximal pvs (i.e., pvs still considered as valid) : "
368                  << mViewCellsManager->GetMaxPvsSize() << endl;
369
370        //-- store rays
371        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
372        {
373                VssRay *ray = *rit;
374
375                float minT, maxT;
376
377                static Ray hray;
378                hray.Init(*ray);
379
380                // TODO: not very efficient to implictly cast between rays types
381                if (mBox.GetRaySegment(hray, minT, maxT))
382                {
383                        float len = ray->Length();
384
385                        if (!len)
386                                len = Limits::Small;
387
388                        rays->push_back(RayInfo(ray, minT / len, maxT / len));
389                }
390        }
391
392        // normalize
393        if (mUseAreaForPvs)
394                mTermMinProbability *= mBox.SurfaceArea();
395        else
396                mTermMinProbability *= mBox.GetVolume();
397
398        // throw out unnecessary polygons
399        PreprocessPolygons(polys);
400
401        mBspStats.polys = (int)polys.size();
402        mGlobalCostMisses = 0;
403
404        cout << "finished" << endl;
405
406        Debug << "\nPolygon extraction: " << (int)polys.size() << " polys extracted from "
407                  << (int)sampleRays.size() << " rays in "
408                  << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;
409
410        // use split cost priority queue
411        if (mUseSplitCostQueue)
412        {
413                ConstructWithSplitQueue(polys, rays);
414        }
415        else
416        {
417                Construct(polys, rays);
418        }
419
420        // clean up polygons
421        CLEAR_CONTAINER(polys);
422}
423
424
425// TODO: return memory usage in MB
426float VspBspTree::GetMemUsage() const
427{
428        return (float)
429                 (sizeof(VspBspTree) +
430                  mBspStats.Leaves() * sizeof(BspLeaf) +
431                  mCreatedViewCells * sizeof(BspViewCell) +
432                  mBspStats.pvs * sizeof(ObjectPvsData) +
433                  mBspStats.Interior() * sizeof(BspInterior) +
434                  mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
435}
436
437
438
439void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
440{
441        VspBspTraversalQueue tQueue;
442
443        mRoot = new BspLeaf();
444
445        // constrruct root node geometry
446        BspNodeGeometry *geom = new BspNodeGeometry();
447        ConstructGeometry(mRoot, *geom);
448
449        const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();
450
451        VspBspTraversalData tData(mRoot,
452                                                          new PolygonContainer(polys),
453                                                          0,
454                                                          rays,
455                              ComputePvsSize(*rays),
456                                                          prop,
457                                                          geom);
458
459        EvalPriority(tData);
460
461
462        // first node is kd node, i.e. an axis aligned box
463        if (1)
464        tData.mIsKdNode = true;
465        else
466                tData.mIsKdNode = false;
467
468        tQueue.push(tData);
469
470
471        mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
472        mTotalPvsSize = tData.mPvs;
473       
474        mSubdivisionStats
475                        << "#ViewCells\n1\n" <<  endl
476                        << "#RenderCostDecrease\n0\n" << endl
477                        << "#TotalRenderCost\n" << mTotalCost << endl
478                        << "#AvgRenderCost\n" << mTotalPvsSize << endl;
479
480        Debug << "total cost: " << mTotalCost << endl;
481       
482       
483        mBspStats.Start();
484        cout << "Constructing vsp bsp tree ... \n";
485
486        long startTime = GetTime();     
487        int nLeaves = 500;
488        int nViewCells = 500;
489
490        // used for intermediate time measurements and progress
491        long interTime = GetTime();     
492
493        mOutOfMemory = false;
494
495        mCreatedViewCells = 0;
496       
497        while (!tQueue.empty())
498        {
499                tData = tQueue.top();
500            tQueue.pop();               
501
502                if (0 && !mOutOfMemory)
503                {
504                        float mem = GetMemUsage();
505
506                        if (mem > mMaxMemory)
507                        {
508                                mOutOfMemory = true;
509                                Debug << "memory limit reached: " << mem << endl;
510                        }
511                }
512
513                // subdivide leaf node
514                BspNode *r = Subdivide(tQueue, tData);
515
516                if (r == mRoot)
517                        Debug << "VSP BSP tree construction time spent at root: "
518                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
519
520                if (mBspStats.Leaves() >= nLeaves)
521                {
522                        nLeaves += 500;
523
524                        cout << "leaves=" << mBspStats.Leaves() << endl;
525                        Debug << "needed "
526                                  << TimeDiff(interTime, GetTime())*1e-3
527                                  << " secs to create 500 view cells" << endl;
528                        interTime = GetTime();
529                }
530
531                if (mCreatedViewCells >= nViewCells)
532                {
533                        nViewCells += 500;
534
535                        cout << "generated " << mCreatedViewCells << " viewcells" << endl;
536                }
537        }
538
539        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
540        cout << "finished\n";
541
542        mBspStats.Stop();
543}
544
545
546
547void VspBspTree::ConstructWithSplitQueue(const PolygonContainer &polys,
548                                                                                          RayInfoContainer *rays)
549{
550        VspBspSplitQueue tQueue;
551
552        mRoot = new BspLeaf();
553
554        // constrruct root node geometry
555        BspNodeGeometry *geom = new BspNodeGeometry();
556        ConstructGeometry(mRoot, *geom);
557
558        const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();
559
560        VspBspTraversalData tData(mRoot,
561                                                          new PolygonContainer(polys),
562                                                          0,
563                                                          rays,
564                              ComputePvsSize(*rays),
565                                                          prop,
566                                                          geom);
567
568
569        // compute first split candidate
570        VspBspSplitCandidate splitCandidate;
571        EvalSplitCandidate(tData, splitCandidate);
572
573        tQueue.push(splitCandidate);
574
575        mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
576        mTotalPvsSize = tData.mPvs;
577       
578        mSubdivisionStats
579                        << "#ViewCells\n1\n" <<  endl
580                        << "#RenderCostDecrease\n0\n" << endl
581                        << "#dummy\n0\n" << endl
582                        << "#TotalRenderCost\n" << mTotalCost << endl
583                        << "#AvgRenderCost\n" << mTotalPvsSize << endl;
584
585        Debug << "total cost: " << mTotalCost << endl;
586       
587       
588        mBspStats.Start();
589        cout << "Constructing vsp bsp tree ... \n";
590
591        long startTime = GetTime();     
592        int nLeaves = 500;
593        int nViewCells = 500;
594
595        // used for intermediate time measurements and progress
596        long interTime = GetTime();     
597
598        mOutOfMemory = false;
599
600        mCreatedViewCells = 0;
601       
602        while (!tQueue.empty())
603        {
604                splitCandidate = tQueue.top();
605            tQueue.pop();               
606
607                // cost ratio of cost decrease / totalCost
608                float costRatio = splitCandidate.GetCost() / mTotalCost;
609
610                //Debug << "cost ratio: " << costRatio << endl;
611
612                if (costRatio < mTermMinGlobalCostRatio)
613                        ++ mGlobalCostMisses;
614               
615                if (0 && !mOutOfMemory)
616                {
617                        float mem = GetMemUsage();
618
619                        if (mem > mMaxMemory)
620                        {
621                                mOutOfMemory = true;
622                                Debug << "memory limit reached: " << mem << endl;
623                        }
624                }
625
626                // subdivide leaf node
627                BspNode *r = Subdivide(tQueue, splitCandidate);
628
629                if (r == mRoot)
630                        Debug << "VSP BSP tree construction time spent at root: "
631                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
632
633                if (mBspStats.Leaves() >= nLeaves)
634                {
635                        nLeaves += 500;
636
637                        cout << "leaves=" << mBspStats.Leaves() << endl;
638                        Debug << "needed "
639                                  << TimeDiff(interTime, GetTime())*1e-3
640                                  << " secs to create 500 view cells" << endl;
641                        interTime = GetTime();
642                }
643
644                if (mCreatedViewCells == nViewCells)
645                {
646                        nViewCells += 500;
647
648                        cout << "generated " << mCreatedViewCells << " viewcells" << endl;
649                }
650        }
651
652        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
653        cout << "finished\n";
654
655        mBspStats.Stop();
656}
657
658
659bool VspBspTree::LocalTerminationCriteriaMet(const VspBspTraversalData &data) const
660{
661        return
662                (((int)data.mRays->size() <= mTermMinRays) ||
663                 (data.mPvs <= mTermMinPvs)   ||
664                 (data.mProbability <= mTermMinProbability) ||
665                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
666                 (data.mDepth >= mTermMaxDepth));
667}
668
669
670bool VspBspTree::GlobalTerminationCriteriaMet(const VspBspTraversalData &data) const
671{
672        return
673                (mOutOfMemory
674                || (mBspStats.Leaves() >= mMaxViewCells)
675                || (mGlobalCostMisses >= mTermGlobalCostMissTolerance)
676                 );
677}
678
679
680
681BspNode *VspBspTree::Subdivide(VspBspTraversalQueue &tQueue,
682                                                           VspBspTraversalData &tData)
683{
684        BspNode *newNode = tData.mNode;
685
686        if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
687        {
688                PolygonContainer coincident;
689
690                VspBspTraversalData tFrontData;
691                VspBspTraversalData tBackData;
692
693                // create new interior node and two leaf nodes
694                // or return leaf as it is (if maxCostRatio missed)
695                int splitAxis;
696                bool splitFurther = true;
697                int maxCostMisses = tData.mMaxCostMisses;
698               
699                Plane3 splitPlane;
700                BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
701               
702                // choose next split plane
703                if (!SelectPlane(splitPlane, leaf, tData, tFrontData, tBackData, splitAxis))
704                {
705                        ++ maxCostMisses;
706
707                        if (maxCostMisses > mTermMissTolerance)
708                        {
709                                // terminate branch because of max cost
710                                ++ mBspStats.maxCostNodes;
711                                splitFurther = false;
712                        }
713                }
714       
715                // if this a valid split => subdivide this node further
716                if (splitFurther) //-- continue subdivision
717                {
718                        newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
719
720                        if (splitAxis < 3)
721                                ++ mBspStats.splits[splitAxis];
722                        else
723                                ++ mBspStats.polySplits;
724
725                        tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && (splitAxis < 3));
726                        tFrontData.mAxis = tBackData.mAxis = splitAxis;
727
728                        // how often was max cost ratio missed in this branch?
729                        tFrontData.mMaxCostMisses = maxCostMisses;
730                        tBackData.mMaxCostMisses = maxCostMisses;
731
732                        EvalPriority(tFrontData);
733                        EvalPriority(tBackData);
734
735                        // evaluate subdivision stats
736                        if (1)
737                        {
738                                float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
739                                float cBack = (float)tBackData.mPvs * tBackData.mProbability;
740                                float cData = (float)tData.mPvs * tData.mProbability;;
741
742                float costDecr = (cFront + cBack - cData) / mBox.GetVolume();
743
744                                mTotalCost += costDecr;
745                                mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
746
747                                mSubdivisionStats
748                                                << "#ViewCells\n" << mBspStats.Leaves() << endl
749                                                << "#RenderCostDecrease\n" << -costDecr << endl
750                                                << "#TotalRenderCost\n" << mTotalCost << endl
751                                                << "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mBspStats.Leaves() << endl;
752                        }
753
754                        // push the children on the stack
755                        tQueue.push(tFrontData);
756                        tQueue.push(tBackData);
757
758                        // delete old leaf node
759                        DEL_PTR(tData.mNode);
760                }
761        }
762
763        //-- terminate traversal and create new view cell
764        if (newNode->IsLeaf())
765        {
766                BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);
767                BspViewCell *viewCell = new BspViewCell();
768               
769                leaf->SetViewCell(viewCell);
770       
771                //-- update pvs
772                int conSamp = 0;
773                float sampCon = 0.0f;
774                AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
775
776                mBspStats.contributingSamples += conSamp;
777                mBspStats.sampleContributions +=(int) sampCon;
778
779                //-- store additional info
780                if (mStoreRays)
781                {
782                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
783                        for (it = tData.mRays->begin(); it != it_end; ++ it)
784                        {
785                                (*it).mRay->Ref();                     
786                                leaf->mVssRays.push_back((*it).mRay);
787                        }
788                }
789
790                // should I check here?
791                if (0 && !mViewCellsManager->CheckValidity(viewCell, 0, mViewCellsManager->GetMaxPvsSize()))
792                {
793                        viewCell->SetValid(false);
794                        leaf->SetTreeValid(false);
795                        PropagateUpValidity(leaf);
796
797                        ++ mBspStats.invalidLeaves;
798                }
799               
800        viewCell->mLeaf = leaf;
801
802                if (mUseAreaForPvs)
803                        viewCell->SetArea(tData.mProbability);
804                else
805                        viewCell->SetVolume(tData.mProbability);
806
807                leaf->mProbability = tData.mProbability;
808
809                EvaluateLeafStats(tData);               
810        }
811
812        //-- cleanup
813        tData.Clear();
814
815        return newNode;
816}
817
818// subdivide using a split plane queue
819BspNode *VspBspTree::Subdivide(VspBspSplitQueue &tQueue,
820                                                           VspBspSplitCandidate &splitCandidate)
821{
822        VspBspTraversalData &tData = splitCandidate.mParentData;
823
824        BspNode *newNode = tData.mNode;
825
826        if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
827        {       
828                PolygonContainer coincident;
829
830                VspBspTraversalData tFrontData;
831                VspBspTraversalData tBackData;
832
833                //-- continue subdivision
834               
835                // create new interior node and two leaf node
836                const Plane3 splitPlane = splitCandidate.mSplitPlane;
837                               
838                newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
839       
840                const int splitAxis = splitCandidate.mSplitAxis;
841                const int maxCostMisses = splitCandidate.mMaxCostMisses;
842
843                if (splitAxis < 3)
844                        ++ mBspStats.splits[splitAxis];
845                else
846                        ++ mBspStats.polySplits;
847
848                tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && (splitAxis < 3));
849                tFrontData.mAxis = tBackData.mAxis = splitAxis;
850
851                // how often was max cost ratio missed in this branch?
852                tFrontData.mMaxCostMisses = maxCostMisses;
853                tBackData.mMaxCostMisses = maxCostMisses;
854                       
855               
856                if (1)
857                {
858                        float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
859                        float cBack = (float)tBackData.mPvs * tBackData.mProbability;
860                        float cData = (float)tData.mPvs * tData.mProbability;
861
862                       
863                        float costDecr =
864                                (cFront + cBack - cData) / mBox.GetVolume();
865
866                        mTotalCost += costDecr;
867                        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
868
869                        mSubdivisionStats
870                                        << "#ViewCells\n" << mBspStats.Leaves() << endl
871                                        << "#RenderCostDecrease\n" << -costDecr << endl
872                                        << "#dummy\n" << splitCandidate.GetCost() << endl
873                                        << "#TotalRenderCost\n" << mTotalCost << endl
874                                        << "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mBspStats.Leaves() << endl;
875                }
876
877       
878                //-- push the new split candidates on the stack
879                VspBspSplitCandidate frontCandidate;
880                VspBspSplitCandidate backCandidate;
881
882                EvalSplitCandidate(tFrontData, frontCandidate);
883                EvalSplitCandidate(tBackData, backCandidate);
884       
885                tQueue.push(frontCandidate);
886                tQueue.push(backCandidate);
887       
888                // delete old leaf node
889                DEL_PTR(tData.mNode);
890        }
891
892
893        //-- terminate traversal and create new view cell
894        if (newNode->IsLeaf())
895        {
896                BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);
897                BspViewCell *viewCell = new BspViewCell();
898
899        leaf->SetViewCell(viewCell);
900               
901                //-- update pvs
902                int conSamp = 0;
903                float sampCon = 0.0f;
904                AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
905
906                mBspStats.contributingSamples += conSamp;
907                mBspStats.sampleContributions +=(int) sampCon;
908
909                //-- store additional info
910                if (mStoreRays)
911                {
912                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
913                        for (it = tData.mRays->begin(); it != it_end; ++ it)
914                        {
915                                (*it).mRay->Ref();                     
916                                leaf->mVssRays.push_back((*it).mRay);
917                        }
918                }
919
920                // should I check here?
921                viewCell->mLeaf = leaf;
922
923                if (mUseAreaForPvs)
924                        viewCell->SetArea(tData.mProbability);
925                else
926                        viewCell->SetVolume(tData.mProbability);
927
928        leaf->mProbability = tData.mProbability;
929
930                EvaluateLeafStats(tData);               
931        }
932
933        //-- cleanup
934        tData.Clear();
935
936        return newNode;
937}
938
939
940void VspBspTree::EvalPriority(VspBspTraversalData &tData) const
941{
942    tData.mPriority = mBreathFirstSplits ?
943                (float)-tData.mDepth : tData.mPvs * tData.mProbability;
944        //cout << "priority: " << tData.mPriority << endl;
945}
946
947
948void VspBspTree::EvalSplitCandidate(VspBspTraversalData &tData,
949                                                                        VspBspSplitCandidate &splitData)
950{
951        VspBspTraversalData frontData;
952        VspBspTraversalData backData;
953
954        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
955
956        // compute locally best split plane
957    bool success = SelectPlane(splitData.mSplitPlane, leaf, tData,
958                                                           frontData, backData, splitData.mSplitAxis);
959
960        // TODO: reuse
961        DEL_PTR(frontData.mGeometry);
962        DEL_PTR(backData.mGeometry);
963       
964        // compute global decrease in render cost
965        splitData.mRenderCost = EvalRenderCostDecrease(splitData.mSplitPlane, tData);
966        splitData.mParentData = tData;
967        splitData.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
968}
969
970
971BspInterior *VspBspTree::SubdivideNode(const Plane3 &splitPlane,
972                                                                           VspBspTraversalData &tData,
973                                                                           VspBspTraversalData &frontData,
974                                                                           VspBspTraversalData &backData,
975                                                                           PolygonContainer &coincident)
976{
977        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
978       
979        //-- the front and back traversal data is filled with the new values
980        frontData.mDepth = tData.mDepth + 1;
981        frontData.mPolygons = new PolygonContainer();
982        frontData.mRays = new RayInfoContainer();
983       
984        backData.mDepth = tData.mDepth + 1;
985        backData.mPolygons = new PolygonContainer();
986        backData.mRays = new RayInfoContainer();
987       
988
989        //-- subdivide rays
990        SplitRays(splitPlane,
991                          *tData.mRays,
992                          *frontData.mRays,
993                          *backData.mRays);
994
995
996        // compute pvs
997        frontData.mPvs = ComputePvsSize(*frontData.mRays);
998        backData.mPvs = ComputePvsSize(*backData.mRays);
999
1000        // split front and back node geometry and compute area
1001       
1002        // if geometry was not already computed
1003        if (!frontData.mGeometry && !backData.mGeometry)
1004        {
1005                frontData.mGeometry = new BspNodeGeometry();
1006                backData.mGeometry = new BspNodeGeometry();
1007
1008                tData.mGeometry->SplitGeometry(*frontData.mGeometry,
1009                                                                           *backData.mGeometry,
1010                                                                           splitPlane,
1011                                                                           mBox,
1012                                                                           //0.0f);
1013                                                                           mEpsilon);
1014               
1015                if (mUseAreaForPvs)
1016                {
1017                        frontData.mProbability = frontData.mGeometry->GetArea();
1018                        backData.mProbability = backData.mGeometry->GetArea();
1019                }
1020                else
1021                {
1022                        frontData.mProbability = frontData.mGeometry->GetVolume();
1023                        backData.mProbability = tData.mProbability - frontData.mProbability;
1024
1025                        if (frontData.mProbability < -0.00001)
1026                                Debug << "fatal error f: " << frontData.mProbability << endl;
1027                        if (backData.mProbability < -0.00001)
1028                                Debug << "fatal error b: " << backData.mProbability << endl;
1029
1030                        // clamp because of precision issues
1031                        if (0)
1032                        {
1033                                if (frontData.mProbability < 0) frontData.mProbability = 0;
1034                                if (backData.mProbability < 0) backData.mProbability = 0;
1035                        }
1036                }
1037        }
1038
1039       
1040    // subdivide polygons
1041        SplitPolygons(splitPlane,
1042                                  *tData.mPolygons,
1043                      *frontData.mPolygons,
1044                                  *backData.mPolygons,
1045                                  coincident);
1046
1047
1048
1049        ///////////////////////////////////////
1050        // subdivide further
1051
1052        // store maximal and minimal depth
1053        if (tData.mDepth > mBspStats.maxDepth)
1054        {
1055                Debug << "max depth increases to " << tData.mDepth << " at " << mBspStats.Leaves() << " leaves" << endl;
1056                mBspStats.maxDepth = tData.mDepth;
1057        }
1058
1059        mBspStats.nodes += 2;
1060
1061   
1062        BspInterior *interior = new BspInterior(splitPlane);
1063
1064#ifdef _DEBUG
1065        Debug << interior << endl;
1066#endif
1067
1068
1069        //-- create front and back leaf
1070
1071        BspInterior *parent = leaf->GetParent();
1072
1073        // replace a link from node's parent
1074        if (parent)
1075        {
1076                parent->ReplaceChildLink(leaf, interior);
1077                interior->SetParent(parent);
1078        }
1079        else // new root
1080        {
1081                mRoot = interior;
1082        }
1083
1084        // and setup child links
1085        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
1086
1087        frontData.mNode = interior->GetFront();
1088        backData.mNode = interior->GetBack();
1089
1090        interior->mTimeStamp = mTimeStamp ++;
1091       
1092
1093        //DEL_PTR(leaf);
1094        return interior;
1095}
1096
1097
1098void VspBspTree::AddToPvs(BspLeaf *leaf,
1099                                                  const RayInfoContainer &rays,
1100                                                  float &sampleContributions,
1101                                                  int &contributingSamples)
1102{
1103  sampleContributions = 0;
1104  contributingSamples = 0;
1105 
1106  RayInfoContainer::const_iterator it, it_end = rays.end();
1107 
1108  ViewCell *vc = leaf->GetViewCell();
1109 
1110  // add contributions from samples to the PVS
1111  for (it = rays.begin(); it != it_end; ++ it)
1112        {
1113          float sc = 0.0f;
1114          VssRay *ray = (*it).mRay;
1115          bool madeContrib = false;
1116          float contribution;
1117          if (ray->mTerminationObject) {
1118                if (vc->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf, contribution))
1119                  madeContrib = true;
1120                sc += contribution;
1121          }
1122         
1123          if (ray->mOriginObject) {
1124                if (vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf, contribution))
1125                  madeContrib = true;
1126                sc += contribution;
1127          }
1128         
1129          sampleContributions += sc;
1130          if (madeContrib)
1131                  ++ contributingSamples;
1132               
1133          //leaf->mVssRays.push_back(ray);
1134        }
1135}
1136
1137
1138void VspBspTree::SortSplitCandidates(const RayInfoContainer &rays,
1139                                                                         const int axis,
1140                                                                         float minBand,
1141                                                                         float maxBand)
1142{
1143        mSplitCandidates->clear();
1144
1145        int requestedSize = 2 * (int)(rays.size());
1146        // creates a sorted split candidates array
1147        if (mSplitCandidates->capacity() > 500000 &&
1148                requestedSize < (int)(mSplitCandidates->capacity() / 10) )
1149        {
1150        delete mSplitCandidates;
1151                mSplitCandidates = new vector<SortableEntry>;
1152        }
1153
1154        mSplitCandidates->reserve(requestedSize);
1155
1156        float pos;
1157
1158        // float values => don't compare with exact values
1159        if (0)
1160        {
1161                minBand += Limits::Small;
1162                maxBand -= Limits::Small;
1163        }
1164
1165        // insert all queries
1166        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
1167        {
1168                const bool positive = (*ri).mRay->HasPosDir(axis);
1169                               
1170                pos = (*ri).ExtrapOrigin(axis);
1171                // clamp to min / max band
1172                if (0) ClipValue(pos, minBand, maxBand);
1173               
1174                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
1175                        pos, (*ri).mRay));
1176
1177                pos = (*ri).ExtrapTermination(axis);
1178                // clamp to min / max band
1179                if (0) ClipValue(pos, minBand, maxBand);
1180
1181                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
1182                        pos, (*ri).mRay));
1183        }
1184
1185        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
1186}
1187
1188
1189float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
1190                                                                                  const AxisAlignedBox3 &box,
1191                                                                                  const int pvsSize,
1192                                                                                  const int axis,
1193                                          float &position)
1194{
1195        const float minBox = box.Min(axis);
1196        const float maxBox = box.Max(axis);
1197        const float sizeBox = maxBox - minBox;
1198
1199        const float minBand = minBox + 0.1f * sizeBox;
1200        const float maxBand = minBox + 0.9f * sizeBox;
1201
1202        SortSplitCandidates(rays, axis, minBand, maxBand);
1203
1204        // go through the lists, count the number of objects left and right
1205        // and evaluate the following cost funcion:
1206        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
1207
1208        int pvsl = 0, pvsr = pvsSize;
1209
1210        int pvsBack = pvsl;
1211        int pvsFront = pvsr;
1212
1213        float sum = (float)pvsSize * sizeBox;
1214        float minSum = 1e20f;
1215
1216        // if no border can be found, take mid split
1217        position = minBox + 0.5f * sizeBox;
1218
1219        Intersectable::NewMail();
1220
1221        RayInfoContainer::const_iterator ri, ri_end = rays.end();
1222
1223        // set all object as belonging to the front pvs
1224        for(ri = rays.begin(); ri != ri_end; ++ ri)
1225        {
1226                Intersectable *oObject = (*ri).mRay->mOriginObject;
1227                Intersectable *tObject = (*ri).mRay->mTerminationObject;
1228
1229                if (oObject)
1230                {
1231                        if (!oObject->Mailed())
1232                        {
1233                                oObject->Mail();
1234                                oObject->mCounter = 1;
1235                        }
1236                        else
1237                        {
1238                                ++ oObject->mCounter;
1239                        }
1240                }
1241                if (tObject)
1242                {
1243                        if (!tObject->Mailed())
1244                        {
1245                                tObject->Mail();
1246                                tObject->mCounter = 1;
1247                        }
1248                        else
1249                        {
1250                                ++ tObject->mCounter;
1251                        }
1252                }
1253        }
1254
1255        Intersectable::NewMail();
1256
1257        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
1258
1259        for (ci = mSplitCandidates->begin(); ci < ci_end; ++ ci)
1260        {
1261                VssRay *ray;
1262                ray = (*ci).ray;
1263               
1264                Intersectable *oObject = ray->mOriginObject;
1265                Intersectable *tObject = ray->mTerminationObject;
1266               
1267
1268                switch ((*ci).type)
1269                {
1270                        case SortableEntry::ERayMin:
1271                                {
1272                                        if (oObject && !oObject->Mailed())
1273                                        {
1274                                                oObject->Mail();
1275                                                ++ pvsl;
1276                                        }
1277                                        if (tObject && !tObject->Mailed())
1278                                        {
1279                                                tObject->Mail();
1280                                                ++ pvsl;
1281                                        }
1282                                        break;
1283                                }
1284                        case SortableEntry::ERayMax:
1285                                {
1286                                        if (oObject)
1287                                        {
1288                                                if (-- oObject->mCounter == 0)
1289                                                        -- pvsr;
1290                                        }
1291
1292                                        if (tObject)
1293                                        {
1294                                                if (-- tObject->mCounter == 0)
1295                                                        -- pvsr;
1296                                        }
1297
1298                                        break;
1299                                }
1300                }
1301
1302                // Note: sufficient to compare size of bounding boxes of front and back side?
1303                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
1304                {
1305                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
1306
1307                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
1308                        // cout<<"cost= "<<sum<<endl;
1309
1310                        if (sum < minSum)
1311                        {
1312                                minSum = sum;
1313                                position = (*ci).value;
1314                               
1315                                pvsBack = pvsl;
1316                                pvsFront = pvsr;
1317                        }
1318                }
1319        }
1320       
1321        //Debug << "pos: " << position << " sizebox: " << sizeBox << " pvs: " << pvsSize << endl;
1322
1323        // -- compute cost
1324        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1325        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1326
1327        const float pOverall = sizeBox;
1328
1329        const float pBack = position - minBox;
1330        const float pFront = maxBox - position;
1331       
1332        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
1333    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1334        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1335       
1336        const float oldRenderCost = penaltyOld * pOverall;
1337        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1338
1339        float ratio = mPvsFactor * newRenderCost / (oldRenderCost + Limits::Small);
1340
1341        //Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1342        //     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
1343
1344        return ratio;
1345}
1346
1347
1348float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
1349                                                                                 const VspBspTraversalData &tData,
1350                                                                                 int &axis,
1351                                                                                 BspNodeGeometry **frontGeom,
1352                                                                                 BspNodeGeometry **backGeom,
1353                                                                                 float &pFront,
1354                                                                                 float &pBack,
1355                                                                                 const bool isKdNode)
1356{
1357        float nPosition[3];
1358        float nCostRatio[3];
1359        float nProbFront[3];
1360        float nProbBack[3];
1361
1362        BspNodeGeometry *nFrontGeom[3];
1363        BspNodeGeometry *nBackGeom[3];
1364
1365        for (int i = 0; i < 3; ++ i)
1366        {
1367                nFrontGeom[i] = NULL;
1368                nBackGeom[i] = NULL;
1369        }
1370
1371        int bestAxis = -1;
1372
1373        // create bounding box of node geometry
1374        AxisAlignedBox3 box;
1375               
1376        //TODO: for kd split geometry already is box => only take minmax vertices
1377        if (1)
1378        {
1379                tData.mGeometry->GetBoundingBox(box);
1380        }
1381        else
1382        {
1383                box.Initialize();
1384                RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();
1385
1386                for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
1387                        box.Include((*ri).ExtrapTermination());
1388        }
1389
1390        int sAxis = 0;
1391
1392        bool useSpecialAxis = mOnlyDrivingAxis || mUseRandomAxis || mSimulateOctree;
1393
1394        // use some kind of specialised fixed axis
1395        if (mOnlyDrivingAxis)
1396                sAxis = box.Size().DrivingAxis();
1397        else if (mUseRandomAxis)
1398                sAxis = Random(3);
1399        else if (mSimulateOctree)
1400                sAxis = (tData.mAxis + 1) % 3;
1401               
1402        //Debug << "use special axis: " << useSpecialAxis << endl;
1403        //Debug << "axis: " << sAxis << " drivingaxis: " << box.Size().DrivingAxis();
1404
1405        for (axis = 0; axis < 3; ++ axis)
1406        {
1407                if (!useSpecialAxis || (axis == sAxis))
1408                {
1409                        if (!mUseCostHeuristics)
1410                        {
1411                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1412                                Vector3 normal(0,0,0); normal[axis] = 1.0f;
1413
1414                                // allows faster split because we have axis aligned kd tree boxes
1415                                if (0 && isKdNode)
1416                                {
1417                                        nCostRatio[axis] = EvalAxisAlignedSplitCost(tData,
1418                                                                                                                                box,
1419                                                                                                                                axis,
1420                                                                                                                                nPosition[axis],
1421                                                                                                                                nProbFront[axis],
1422                                                                                                                                nProbBack[axis]);
1423                                       
1424                                        Vector3 pos;
1425                                       
1426                                        // create back geometry from box
1427                                        pos = box.Max(); pos[axis] = nPosition[axis];
1428                                        AxisAlignedBox3 bBox(box.Min(), pos);
1429                                        PolygonContainer fPolys;
1430                                        bBox.ExtractPolys(fPolys);
1431
1432                                        nBackGeom[axis] = new BspNodeGeometry(fPolys);
1433       
1434                                        // create front geometry from box
1435                                        pos = box.Min(); pos[axis] = nPosition[axis];
1436                                        AxisAlignedBox3 fBox(pos, box.Max());
1437
1438                                        PolygonContainer bPolys;
1439                                        fBox.ExtractPolys(bPolys);
1440                                        nFrontGeom[axis] = new BspNodeGeometry(bPolys);
1441                                }
1442                                else
1443                                {
1444                                        nFrontGeom[axis] = new BspNodeGeometry();
1445                                        nBackGeom[axis] = new BspNodeGeometry();
1446
1447                                        nCostRatio[axis] =
1448                                                EvalSplitPlaneCost(Plane3(normal, nPosition[axis]),
1449                                                                                   tData, *nFrontGeom[axis], *nBackGeom[axis],
1450                                                                                   nProbFront[axis], nProbBack[axis]);
1451                                }
1452                        }
1453                        else
1454                        {
1455                                nCostRatio[axis] =
1456                                        BestCostRatioHeuristics(*tData.mRays,
1457                                                                                    box,
1458                                                                                        tData.mPvs,
1459                                                                                        axis,
1460                                                                                        nPosition[axis]);
1461                        }
1462
1463                        if (bestAxis == -1)
1464                        {
1465                                bestAxis = axis;
1466                        }
1467                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1468                        {
1469                                bestAxis = axis;
1470                        }
1471
1472                }
1473        }
1474
1475        //-- assign values
1476        axis = bestAxis;
1477        pFront = nProbFront[bestAxis];
1478        pBack = nProbBack[bestAxis];
1479
1480        // assign best split nodes geometry
1481        *frontGeom = nFrontGeom[bestAxis];
1482        *backGeom = nBackGeom[bestAxis];
1483
1484        // and delete other geometry
1485        DEL_PTR(nFrontGeom[(bestAxis + 1) % 3]);
1486        DEL_PTR(nBackGeom[(bestAxis + 2) % 3]);
1487
1488        //-- split plane
1489    Vector3 normal(0,0,0); normal[bestAxis] = 1;
1490        plane = Plane3(normal, nPosition[bestAxis]);
1491
1492        return nCostRatio[bestAxis];
1493}
1494
1495
1496bool VspBspTree::SelectPlane(Plane3 &bestPlane,
1497                                                         BspLeaf *leaf,
1498                                                         VspBspTraversalData &data,                                                     
1499                                                         VspBspTraversalData &frontData,
1500                                                         VspBspTraversalData &backData,
1501                                                         int &splitAxis)
1502{
1503        // simplest strategy: just take next polygon
1504        if (mSplitPlaneStrategy & RANDOM_POLYGON)
1505        {
1506        if (!data.mPolygons->empty())
1507                {
1508                        const int randIdx =
1509                                (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
1510                        Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
1511
1512                        bestPlane = nextPoly->GetSupportingPlane();
1513                        return true;
1514                }
1515        }
1516
1517        //-- use heuristics to find appropriate plane
1518
1519        // intermediate plane
1520        Plane3 plane;
1521        float lowestCost = MAX_FLOAT;
1522       
1523        // decides if the first few splits should be only axisAligned
1524        const bool onlyAxisAligned  =
1525                (mSplitPlaneStrategy & AXIS_ALIGNED) &&
1526                ((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
1527                ((int)data.GetAvgRayContribution() < mTermMaxRayContriForAxisAligned);
1528       
1529        const int limit = onlyAxisAligned ? 0 :
1530                Min((int)data.mPolygons->size(), mMaxPolyCandidates);
1531
1532        float candidateCost;
1533
1534        int maxIdx = (int)data.mPolygons->size();
1535
1536        for (int i = 0; i < limit; ++ i)
1537        {
1538                // the already taken candidates are stored behind maxIdx
1539                // => assure that no index is taken twice
1540                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
1541                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
1542
1543                // swap candidate to the end to avoid testing same plane
1544                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
1545                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
1546
1547                // evaluate current candidate
1548                BspNodeGeometry fGeom, bGeom;
1549                float fArea, bArea;
1550                plane = poly->GetSupportingPlane();
1551                candidateCost = EvalSplitPlaneCost(plane, data, fGeom, bGeom, fArea, bArea);
1552               
1553                if (candidateCost < lowestCost)
1554                {
1555                        bestPlane = plane;
1556                        lowestCost = candidateCost;
1557                }
1558        }
1559
1560        //-- evaluate axis aligned splits
1561        int axis;
1562        BspNodeGeometry *fGeom, *bGeom;
1563        float pFront, pBack;
1564
1565        candidateCost = 99999999.0f;
1566
1567        // option: axis aligned split only if no polygon available
1568        if (!mUsePolygonSplitIfAvailable || data.mPolygons->empty())
1569        {
1570                candidateCost = SelectAxisAlignedPlane(plane,
1571                                                                                           data,
1572                                                                                           axis,
1573                                                                                           &fGeom,
1574                                                                                           &bGeom,
1575                                                                                           pFront,
1576                                                                                           pBack,
1577                                                                                           data.mIsKdNode);     
1578        }
1579
1580        splitAxis = 3;
1581
1582        if (candidateCost < lowestCost)
1583        {       
1584                bestPlane = plane;
1585                lowestCost = candidateCost;
1586                splitAxis = axis;
1587       
1588                // assign already computed values
1589                // we can do this because we always save the
1590                // computed values from the axis aligned splits         
1591
1592                if (fGeom && bGeom)
1593                {
1594                        frontData.mGeometry = fGeom;
1595                        backData.mGeometry = bGeom;
1596       
1597                        frontData.mProbability = pFront;
1598                        backData.mProbability = pBack;
1599                }
1600        }
1601        else
1602        {
1603                DEL_PTR(fGeom);
1604                DEL_PTR(bGeom);
1605        }
1606   
1607
1608//      if (lowestCost > 10)
1609//              Debug << "warning!! lowest cost: " << lowestCost << endl;
1610       
1611#ifdef _DEBUG
1612        Debug << "plane lowest cost: " << lowestCost << endl;
1613#endif
1614
1615        if (lowestCost > mTermMaxCostRatio)
1616        {
1617                return false;
1618        }
1619
1620        return true;
1621}
1622
1623
1624Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
1625{
1626        const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1627
1628        const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
1629        const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();
1630
1631        const Vector3 pt = (maxPt + minPt) * 0.5;
1632        const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());
1633
1634        return Plane3(normal, pt);
1635}
1636
1637
1638Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
1639{
1640        Vector3 pt[3];
1641
1642        int idx[3];
1643        int cmaxT = 0;
1644        int cminT = 0;
1645        bool chooseMin = false;
1646
1647        for (int j = 0; j < 3; ++ j)
1648        {
1649                idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));
1650
1651                if (idx[j] >= (int)rays.size())
1652                {
1653                        idx[j] -= (int)rays.size();
1654
1655                        chooseMin = (cminT < 2);
1656                }
1657                else
1658                        chooseMin = (cmaxT < 2);
1659
1660                RayInfo rayInf = rays[idx[j]];
1661                pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
1662        }
1663
1664        return Plane3(pt[0], pt[1], pt[2]);
1665}
1666
1667
1668Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
1669{
1670        Vector3 pt[3];
1671
1672        int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1673        int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1674
1675        // check if rays different
1676        if (idx1 == idx2)
1677                idx2 = (idx2 + 1) % (int)rays.size();
1678
1679        const RayInfo ray1 = rays[idx1];
1680        const RayInfo ray2 = rays[idx2];
1681
1682        // normal vector of the plane parallel to both lines
1683        const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));
1684
1685        // vector from line 1 to line 2
1686        const Vector3 vd = ray2.ExtrapOrigin() - ray1.ExtrapOrigin();
1687
1688        // project vector on normal to get distance
1689        const float dist = DotProd(vd, norm);
1690
1691        // point on plane lies halfway between the two planes
1692        const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;
1693
1694        return Plane3(norm, planePt);
1695}
1696
1697
1698inline void VspBspTree::GenerateUniqueIdsForPvs()
1699{
1700        Intersectable::NewMail(); sBackId = Intersectable::sMailId;
1701        Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
1702        Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
1703}
1704
1705
1706float VspBspTree::EvalRenderCostDecrease(const Plane3 &candidatePlane,
1707                                                                                 const VspBspTraversalData &data) const
1708{
1709        float pvsFront = 0;
1710        float pvsBack = 0;
1711        float totalPvs = 0;
1712
1713        // probability that view point lies in back / front node
1714        float pOverall = data.mProbability;
1715        float pFront = 0;
1716        float pBack = 0;
1717
1718
1719        // create unique ids for pvs heuristics
1720        GenerateUniqueIdsForPvs();
1721       
1722        for (int i = 0; i < data.mRays->size(); ++ i)
1723        {
1724                RayInfo rayInf = (*data.mRays)[i];
1725
1726                float t;
1727                VssRay *ray = rayInf.mRay;
1728                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
1729
1730                // find front and back pvs for origing and termination object
1731                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1732                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1733        }
1734
1735
1736        BspNodeGeometry geomFront;
1737        BspNodeGeometry geomBack;
1738
1739        // construct child geometry with regard to the candidate split plane
1740        data.mGeometry->SplitGeometry(geomFront,
1741                                                                  geomBack,
1742                                                                  candidatePlane,
1743                                                                  mBox,
1744                                                                  //0.0f);
1745                                                                  mEpsilon);
1746
1747        if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
1748        {
1749                pFront = geomFront.GetVolume();
1750                pBack = pOverall - pFront;
1751
1752                // something is wrong with the volume
1753                if ((pFront < 0.0) || (pBack < 0.0))
1754                {
1755                        Debug << "vol f " << pFront << " b " << pBack << " p " << pOverall << endl;
1756                        Debug << "real vol f " << pFront << " b " << geomBack.GetVolume() << " p " << pOverall << endl;
1757                        Debug << "polys f " << geomFront.Size() << " b " << geomBack.Size() << " data " << data.mGeometry->Size() << endl;
1758                }
1759
1760                // clamp because of possible precision issues
1761                if (0)
1762                {
1763                        if (pFront < 0) pFront = 0;
1764                        if (pBack < 0) pBack = 0;
1765                }
1766        }
1767        else
1768        {
1769                pFront = geomFront.GetArea();
1770                pBack = geomBack.GetArea();
1771        }
1772       
1773
1774        // -- pvs rendering heuristics
1775        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1776        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1777
1778        // only render cost heuristics or combined with standard deviation
1779        const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
1780    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1781        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1782                       
1783        const float oldRenderCost = pOverall * penaltyOld;
1784        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1785
1786        //Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
1787        return (oldRenderCost - newRenderCost) / mBox.GetVolume();
1788}
1789
1790
1791float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
1792                                                                         const VspBspTraversalData &data,
1793                                                                         BspNodeGeometry &geomFront,
1794                                                                         BspNodeGeometry &geomBack,
1795                                                                         float &pFront,
1796                                                                         float &pBack) const
1797{
1798        float cost = 0;
1799
1800        float totalPvs = 0;
1801        float pvsFront = 0;
1802        float pvsBack = 0;
1803       
1804        // probability that view point lies in back / front node
1805        float pOverall = 0;
1806        pFront = 0;
1807        pBack = 0;
1808
1809
1810        int limit;
1811        bool useRand;
1812
1813        // create unique ids for pvs heuristics
1814        GenerateUniqueIdsForPvs();
1815
1816        // choose test rays randomly if too much
1817        if ((int)data.mRays->size() > mMaxTests)
1818        {
1819                useRand = true;
1820                limit = mMaxTests;
1821        }
1822        else
1823        {
1824                useRand = false;
1825                limit = (int)data.mRays->size();
1826        }
1827       
1828        for (int i = 0; i < limit; ++ i)
1829        {
1830                const int testIdx = useRand ?
1831                        (int)RandomValue(0, (Real)((int)data.mRays->size() - 1)) : i;
1832                RayInfo rayInf = (*data.mRays)[testIdx];
1833
1834                float t;
1835                VssRay *ray = rayInf.mRay;
1836                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
1837
1838                // find front and back pvs for origing and termination object
1839                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1840                AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1841        }
1842
1843        // construct child geometry with regard to the candidate split plane
1844        bool splitSuccessFull = data.mGeometry->SplitGeometry(geomFront,
1845                                                                                                                  geomBack,
1846                                                                                                                  candidatePlane,
1847                                                                                                                  mBox,
1848                                                                                                                  //0.0f);
1849                                                                                                                  mEpsilon);
1850
1851        pOverall = data.mProbability;
1852
1853        if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
1854        {
1855                pFront = geomFront.GetVolume();
1856                pBack = pOverall - pFront;
1857               
1858                // HACK: precision issues possible for unbalanced split => don't take this split!
1859                if (1 &&
1860                        (!splitSuccessFull || (pFront <= 0) || (pBack <= 0) ||
1861                        !geomFront.Valid() || !geomBack.Valid()))
1862                {
1863                        Debug << "error f: " << pFront << " b: " << pBack << endl;
1864                        return 99999.9f;
1865                }
1866        }
1867        else
1868        {
1869                pFront = geomFront.GetArea();
1870                pBack = geomBack.GetArea();
1871        }
1872       
1873
1874        // -- pvs rendering heuristics
1875        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1876        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1877
1878        // only render cost heuristics or combined with standard deviation
1879        const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
1880    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1881        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1882                       
1883        const float oldRenderCost = pOverall * penaltyOld;
1884        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1885
1886        float oldCost, newCost;
1887
1888        // only render cost
1889        if (1)
1890        {
1891                oldCost = oldRenderCost;
1892                newCost = newRenderCost;
1893        }
1894        else // also considering standard deviation
1895        {
1896                // standard deviation is difference of back and front pvs
1897                const float expectedCost = 0.5f * (penaltyFront + penaltyBack);
1898
1899                const float newDeviation = 0.5f *
1900                        fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);
1901
1902                const float oldDeviation = penaltyOld;
1903
1904                newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
1905                oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
1906        }
1907
1908        cost += mPvsFactor * newCost / (oldCost + Limits::Small);
1909               
1910
1911#ifdef _DEBUG
1912        Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
1913                  << " frontpvs: " << pvsFront << " pFront: " << pFront
1914                  << " backpvs: " << pvsBack << " pBack: " << pBack << endl << endl;
1915        Debug << "cost: " << cost << endl;
1916#endif
1917
1918        return cost;
1919}
1920
1921
1922int VspBspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box,
1923                                                                                ViewCellContainer &viewCells) const
1924{
1925        stack<bspNodePair> nodeStack;
1926        BspNodeGeometry *rgeom = new BspNodeGeometry();
1927
1928        ConstructGeometry(mRoot, *rgeom);
1929
1930        nodeStack.push(bspNodePair(mRoot, rgeom));
1931 
1932        ViewCell::NewMail();
1933
1934        while (!nodeStack.empty())
1935        {
1936                BspNode *node = nodeStack.top().first;
1937                BspNodeGeometry *geom = nodeStack.top().second;
1938                nodeStack.pop();
1939
1940                const int side = geom->ComputeIntersection(box);
1941               
1942                switch (side)
1943                {
1944                case -1:
1945                        // node geometry is contained in box
1946                        CollectViewCells(node, true, viewCells, true);
1947                        break;
1948
1949                case 0:
1950                        if (node->IsLeaf())
1951                        {
1952                                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1953                       
1954                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
1955                                {
1956                                        leaf->GetViewCell()->Mail();
1957                                        viewCells.push_back(leaf->GetViewCell());
1958                                }
1959                        }
1960                        else
1961                        {
1962                                BspInterior *interior = dynamic_cast<BspInterior *>(node);
1963                       
1964                                BspNode *first = interior->GetFront();
1965                                BspNode *second = interior->GetBack();
1966           
1967                                BspNodeGeometry *firstGeom = new BspNodeGeometry();
1968                                BspNodeGeometry *secondGeom = new BspNodeGeometry();
1969
1970                                geom->SplitGeometry(*firstGeom,
1971                                                                        *secondGeom,
1972                                                                        interior->GetPlane(),
1973                                                                        mBox,
1974                                                                        //0.0000001f);
1975                                                                        mEpsilon);
1976
1977                                nodeStack.push(bspNodePair(first, firstGeom));
1978                                nodeStack.push(bspNodePair(second, secondGeom));
1979                        }
1980                       
1981                        break;
1982                default:
1983                        // default: cull
1984                        break;
1985                }
1986               
1987                DEL_PTR(geom);
1988               
1989        }
1990
1991        return (int)viewCells.size();
1992}
1993
1994
1995float VspBspTree::EvalAxisAlignedSplitCost(const VspBspTraversalData &data,
1996                                                                                   const AxisAlignedBox3 &box,
1997                                                                                   const int axis,
1998                                                                                   const float &position,                                                                                 
1999                                                                                   float &pFront,
2000                                                                                   float &pBack) const
2001{
2002        float pvsTotal = 0;
2003        float pvsFront = 0;
2004        float pvsBack = 0;
2005       
2006        // create unique ids for pvs heuristics
2007        GenerateUniqueIdsForPvs();
2008
2009        const int pvsSize = data.mPvs;
2010
2011        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
2012
2013        // this is the main ray classification loop!
2014        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
2015        {
2016                //if (!(*rit).mRay->IsActive()) continue;
2017
2018                // determine the side of this ray with respect to the plane
2019                float t;
2020                const int side = (*rit).ComputeRayIntersection(axis, position, t);
2021       
2022                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
2023                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
2024        }
2025
2026        //-- pvs heuristics
2027        float pOverall;
2028
2029        //-- compute heurstics
2030        //   we take simplified computation for mid split
2031               
2032        pOverall = data.mProbability;
2033
2034        if (!mUseAreaForPvs)
2035        {   // volume
2036                pBack = pFront = pOverall * 0.5f;
2037               
2038#if 0
2039                // box length substitute for probability
2040                const float minBox = box.Min(axis);
2041                const float maxBox = box.Max(axis);
2042
2043                pBack = position - minBox;
2044                pFront = maxBox - position;
2045                pOverall = maxBox - minBox;
2046#endif
2047        }
2048        else //-- area substitute for probability
2049        {
2050                const int axis2 = (axis + 1) % 3;
2051                const int axis3 = (axis + 2) % 3;
2052
2053                const float faceArea =
2054                        (box.Max(axis2) - box.Min(axis2)) *
2055                        (box.Max(axis3) - box.Min(axis3));
2056
2057                pBack = pFront = pOverall * 0.5f + faceArea;
2058        }
2059
2060#ifdef _DEBUG
2061        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
2062        Debug << pFront << " " << pBack << " " << pOverall << endl;
2063#endif
2064
2065       
2066        const float newCost = pvsBack * pBack + pvsFront * pFront;
2067        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
2068
2069        return  (mCtDivCi + newCost) / oldCost;
2070}
2071
2072
2073void VspBspTree::AddObjToPvs(Intersectable *obj,
2074                                                         const int cf,
2075                                                         float &frontPvs,
2076                                                         float &backPvs,
2077                                                         float &totalPvs) const
2078{
2079        if (!obj)
2080                return;
2081
2082        // new object
2083        if ((obj->mMailbox != sFrontId) &&
2084                (obj->mMailbox != sBackId) &&
2085                (obj->mMailbox != sFrontAndBackId))
2086        {
2087                totalPvs += mViewCellsManager->EvalRenderCost(obj);
2088        }
2089
2090        // TODO: does this really belong to no pvs?
2091        //if (cf == Ray::COINCIDENT) return;
2092
2093        // object belongs to both PVS
2094        if (cf >= 0)
2095        {
2096                if ((obj->mMailbox != sFrontId) &&
2097                        (obj->mMailbox != sFrontAndBackId))
2098                {
2099                        frontPvs += mViewCellsManager->EvalRenderCost(obj);
2100               
2101                        if (obj->mMailbox == sBackId)
2102                                obj->mMailbox = sFrontAndBackId;
2103                        else
2104                                obj->mMailbox = sFrontId;
2105                }
2106        }
2107
2108        if (cf <= 0)
2109        {
2110                if ((obj->mMailbox != sBackId) &&
2111                        (obj->mMailbox != sFrontAndBackId))
2112                {
2113                        backPvs += mViewCellsManager->EvalRenderCost(obj);;
2114               
2115                        if (obj->mMailbox == sFrontId)
2116                                obj->mMailbox = sFrontAndBackId;
2117                        else
2118                                obj->mMailbox = sBackId;
2119                }
2120        }
2121}
2122
2123
2124void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves,
2125                                                           const bool onlyUnmailed,
2126                                                           const int maxPvsSize) const
2127{
2128        stack<BspNode *> nodeStack;
2129        nodeStack.push(mRoot);
2130
2131        while (!nodeStack.empty())
2132        {
2133                BspNode *node = nodeStack.top();
2134                nodeStack.pop();
2135               
2136                if (node->IsLeaf())
2137                {
2138                        // test if this leaf is in valid view space
2139                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2140                        if (leaf->TreeValid() &&
2141                                (!onlyUnmailed || !leaf->Mailed()) &&
2142                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
2143                        {
2144                                leaves.push_back(leaf);
2145                        }
2146                }
2147                else
2148                {
2149                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2150
2151                        nodeStack.push(interior->GetBack());
2152                        nodeStack.push(interior->GetFront());
2153                }
2154        }
2155}
2156
2157
2158AxisAlignedBox3 VspBspTree::GetBoundingBox() const
2159{
2160        return mBox;
2161}
2162
2163
2164BspNode *VspBspTree::GetRoot() const
2165{
2166        return mRoot;
2167}
2168
2169
2170void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
2171{
2172        // the node became a leaf -> evaluate stats for leafs
2173        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
2174
2175
2176        if (data.mPvs > mBspStats.maxPvs)
2177        {
2178                mBspStats.maxPvs = data.mPvs;
2179        }
2180
2181        mBspStats.pvs += data.mPvs;
2182
2183        if (data.mDepth < mBspStats.minDepth)
2184        {
2185                mBspStats.minDepth = data.mDepth;
2186        }
2187       
2188        if (data.mDepth >= mTermMaxDepth)
2189        {
2190                ++ mBspStats.maxDepthNodes;
2191        }
2192
2193        // accumulate rays to compute rays /  leaf
2194        mBspStats.accumRays += (int)data.mRays->size();
2195
2196        if (data.mPvs < mTermMinPvs)
2197                ++ mBspStats.minPvsNodes;
2198
2199        if ((int)data.mRays->size() < mTermMinRays)
2200                ++ mBspStats.minRaysNodes;
2201
2202        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
2203                ++ mBspStats.maxRayContribNodes;
2204
2205        if (data.mProbability <= mTermMinProbability)
2206                ++ mBspStats.minProbabilityNodes;
2207       
2208        // accumulate depth to compute average depth
2209        mBspStats.accumDepth += data.mDepth;
2210
2211        ++ mCreatedViewCells;
2212
2213#ifdef _DEBUG
2214        Debug << "BSP stats: "
2215                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
2216                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
2217          //              << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
2218                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
2219                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
2220                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
2221#endif
2222}
2223
2224
2225int VspBspTree::CastRay(Ray &ray)
2226{
2227        int hits = 0;
2228
2229        stack<BspRayTraversalData> tQueue;
2230
2231        float maxt, mint;
2232
2233        if (!mBox.GetRaySegment(ray, mint, maxt))
2234                return 0;
2235
2236        Intersectable::NewMail();
2237        ViewCell::NewMail();
2238        Vector3 entp = ray.Extrap(mint);
2239        Vector3 extp = ray.Extrap(maxt);
2240
2241        BspNode *node = mRoot;
2242        BspNode *farChild = NULL;
2243
2244        while (1)
2245        {
2246                if (!node->IsLeaf())
2247                {
2248                        BspInterior *in = dynamic_cast<BspInterior *>(node);
2249
2250                        Plane3 splitPlane = in->GetPlane();
2251                        const int entSide = splitPlane.Side(entp);
2252                        const int extSide = splitPlane.Side(extp);
2253
2254                        if (entSide < 0)
2255                        {
2256                                node = in->GetBack();
2257
2258                                if(extSide <= 0) // plane does not split ray => no far child
2259                                        continue;
2260
2261                                farChild = in->GetFront(); // plane splits ray
2262
2263                        } else if (entSide > 0)
2264                        {
2265                                node = in->GetFront();
2266
2267                                if (extSide >= 0) // plane does not split ray => no far child
2268                                        continue;
2269
2270                                farChild = in->GetBack(); // plane splits ray
2271                        }
2272                        else // ray and plane are coincident
2273                        {
2274                                // WHAT TO DO IN THIS CASE ?
2275                                //break;
2276                                node = in->GetFront();
2277                                continue;
2278                        }
2279
2280                        // push data for far child
2281                        tQueue.push(BspRayTraversalData(farChild, extp, maxt));
2282
2283                        // find intersection of ray segment with plane
2284                        float t;
2285                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
2286                        maxt *= t;
2287
2288                } else // reached leaf => intersection with view cell
2289                {
2290                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2291
2292                        if (!leaf->GetViewCell()->Mailed())
2293                        {
2294                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
2295                                leaf->GetViewCell()->Mail();
2296                                ++ hits;
2297                        }
2298
2299                        //-- fetch the next far child from the stack
2300                        if (tQueue.empty())
2301                                break;
2302
2303                        entp = extp;
2304                        mint = maxt; // NOTE: need this?
2305
2306                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
2307                                break;
2308
2309                        BspRayTraversalData &s = tQueue.top();
2310
2311                        node = s.mNode;
2312                        extp = s.mExitPoint;
2313                        maxt = s.mMaxT;
2314
2315                        tQueue.pop();
2316                }
2317        }
2318
2319        return hits;
2320}
2321
2322
2323void VspBspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
2324{
2325        ViewCell::NewMail();
2326       
2327        CollectViewCells(mRoot, onlyValid, viewCells, true);
2328}
2329
2330
2331void VspBspTree::CollapseViewCells()
2332{
2333// TODO
2334#if HAS_TO_BE_REDONE
2335        stack<BspNode *> nodeStack;
2336
2337        if (!mRoot)
2338                return;
2339
2340        nodeStack.push(mRoot);
2341       
2342        while (!nodeStack.empty())
2343        {
2344                BspNode *node = nodeStack.top();
2345                nodeStack.pop();
2346               
2347                if (node->IsLeaf())
2348        {
2349                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
2350
2351                        if (!viewCell->GetValid())
2352                        {
2353                                BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
2354       
2355                                ViewCellContainer leaves;
2356                                mViewCellsTree->CollectLeaves(viewCell, leaves);
2357
2358                                ViewCellContainer::const_iterator it, it_end = leaves.end();
2359
2360                                for (it = leaves.begin(); it != it_end; ++ it)
2361                                {
2362                                        BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
2363                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
2364                                        ++ mBspStats.invalidLeaves;
2365                                }
2366
2367                                // add to unbounded view cell
2368                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
2369                                DEL_PTR(viewCell);
2370                        }
2371                }
2372                else
2373                {
2374                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2375               
2376                        nodeStack.push(interior->GetFront());
2377                        nodeStack.push(interior->GetBack());
2378                }
2379        }
2380
2381        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
2382#endif
2383}
2384
2385
2386void VspBspTree::CollectRays(VssRayContainer &rays)
2387{
2388        vector<BspLeaf *> leaves;
2389
2390        vector<BspLeaf *>::const_iterator lit, lit_end = leaves.end();
2391
2392        for (lit = leaves.begin(); lit != lit_end; ++ lit)
2393        {
2394                BspLeaf *leaf = *lit;
2395                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
2396
2397                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
2398                        rays.push_back(*rit);
2399        }
2400}
2401
2402
2403void VspBspTree::ValidateTree()
2404{
2405        stack<BspNode *> nodeStack;
2406
2407        if (!mRoot)
2408                return;
2409
2410        nodeStack.push(mRoot);
2411       
2412        mBspStats.invalidLeaves = 0;
2413        while (!nodeStack.empty())
2414        {
2415                BspNode *node = nodeStack.top();
2416                nodeStack.pop();
2417               
2418                if (node->IsLeaf())
2419                {
2420                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2421
2422                        if (!leaf->GetViewCell()->GetValid())
2423                                ++ mBspStats.invalidLeaves;
2424
2425                        // validity flags don't match => repair
2426                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
2427                        {
2428                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
2429                                PropagateUpValidity(leaf);
2430                        }
2431                }
2432                else
2433                {
2434                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2435               
2436                        nodeStack.push(interior->GetFront());
2437                        nodeStack.push(interior->GetBack());
2438                }
2439        }
2440
2441        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
2442}
2443
2444
2445
2446void VspBspTree::CollectViewCells(BspNode *root,
2447                                                                  bool onlyValid,
2448                                                                  ViewCellContainer &viewCells,
2449                                                                  bool onlyUnmailed) const
2450{
2451        stack<BspNode *> nodeStack;
2452
2453        if (!root)
2454                return;
2455
2456        nodeStack.push(root);
2457       
2458        while (!nodeStack.empty())
2459        {
2460                BspNode *node = nodeStack.top();
2461                nodeStack.pop();
2462               
2463                if (node->IsLeaf())
2464                {
2465                        if (!onlyValid || node->TreeValid())
2466                        {
2467                                ViewCell *leafVc = dynamic_cast<BspLeaf *>(node)->GetViewCell();
2468
2469                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
2470                                               
2471                                if (!onlyUnmailed || !viewCell->Mailed())
2472                                {
2473                                        viewCell->Mail();
2474                                        viewCells.push_back(viewCell);
2475                                }
2476                        }
2477                }
2478                else
2479                {
2480                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2481               
2482                        nodeStack.push(interior->GetFront());
2483                        nodeStack.push(interior->GetBack());
2484                }
2485        }
2486
2487}
2488
2489
2490void VspBspTree::PreprocessPolygons(PolygonContainer &polys)
2491{
2492        // preprocess: throw out polygons coincident to the view space box (not needed)
2493        PolygonContainer boxPolys;
2494        mBox.ExtractPolys(boxPolys);
2495        vector<Plane3> boxPlanes;
2496
2497        PolygonContainer::iterator pit, pit_end = boxPolys.end();
2498
2499        // extract planes of box
2500        // TODO: can be done more elegantly than first extracting polygons
2501        // and take their planes
2502        for (pit = boxPolys.begin(); pit != pit_end; ++ pit)
2503        {
2504                boxPlanes.push_back((*pit)->GetSupportingPlane());
2505        }
2506
2507        pit_end = polys.end();
2508
2509        for (pit = polys.begin(); pit != pit_end; ++ pit)
2510        {
2511                vector<Plane3>::const_iterator bit, bit_end = boxPlanes.end();
2512               
2513                for (bit = boxPlanes.begin(); (bit != bit_end) && (*pit); ++ bit)
2514                {
2515                        const int cf = (*pit)->ClassifyPlane(*bit, mEpsilon);
2516
2517                        if (cf == Polygon3::COINCIDENT)
2518                        {
2519                                DEL_PTR(*pit);
2520                                //Debug << "coincident!!" << endl;
2521                        }
2522                }
2523        }
2524
2525        // remove deleted entries
2526        for (int i = 0; i < (int)polys.size(); ++ i)
2527        {
2528                while (!polys[i] && (i < (int)polys.size()))
2529                {
2530                        swap(polys[i], polys.back());
2531                        polys.pop_back();
2532                }
2533        }
2534}
2535
2536
2537float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
2538{
2539        float len = 0;
2540
2541        RayInfoContainer::const_iterator it, it_end = rays.end();
2542
2543        for (it = rays.begin(); it != it_end; ++ it)
2544                len += (*it).SegmentLength();
2545
2546        return len;
2547}
2548
2549
2550int VspBspTree::SplitRays(const Plane3 &plane,
2551                                                  RayInfoContainer &rays,
2552                                                  RayInfoContainer &frontRays,
2553                                                  RayInfoContainer &backRays) const
2554{
2555        int splits = 0;
2556
2557        RayInfoContainer::const_iterator it, it_end = rays.end();
2558
2559        for (it = rays.begin(); it != it_end; ++ it)
2560        {
2561                RayInfo bRay = *it;
2562               
2563                VssRay *ray = bRay.mRay;
2564                float t;
2565
2566                // get classification and receive new t
2567                const int cf = bRay.ComputeRayIntersection(plane, t);
2568
2569                switch (cf)
2570                {
2571                case -1:
2572                        backRays.push_back(bRay);
2573                        break;
2574                case 1:
2575                        frontRays.push_back(bRay);
2576                        break;
2577                case 0:
2578                        {
2579                                //-- split ray
2580                                //   test if start point behind or in front of plane
2581                                const int side = plane.Side(bRay.ExtrapOrigin());
2582
2583                                if (side <= 0)
2584                                {
2585                                        backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2586                                        frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2587                                }
2588                                else
2589                                {
2590                                        frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2591                                        backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2592                                }
2593                        }
2594                        break;
2595                default:
2596                        Debug << "Should not come here" << endl;
2597                        break;
2598                }
2599        }
2600
2601        return splits;
2602}
2603
2604
2605void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
2606{
2607        BspNode *lastNode;
2608
2609        do
2610        {
2611                lastNode = n;
2612
2613                // want to get planes defining geometry of this node => don't take
2614                // split plane of node itself
2615                n = n->GetParent();
2616
2617                if (n)
2618                {
2619                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
2620                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
2621
2622            if (interior->GetBack() != lastNode)
2623                                halfSpace.ReverseOrientation();
2624
2625                        halfSpaces.push_back(halfSpace);
2626                }
2627        }
2628        while (n);
2629}
2630
2631
2632void VspBspTree::ConstructGeometry(BspNode *n,
2633                                                                   BspNodeGeometry &geom) const
2634{
2635        vector<Plane3> halfSpaces;
2636        ExtractHalfSpaces(n, halfSpaces);
2637
2638        PolygonContainer candidatePolys;
2639        vector<Plane3> candidatePlanes;
2640
2641        vector<Plane3>::const_iterator pit, pit_end = halfSpaces.end();
2642
2643        // bounded planes are added to the polygons
2644        for (pit = halfSpaces.begin(); pit != pit_end; ++ pit)
2645        {
2646                Polygon3 *p = GetBoundingBox().CrossSection(*pit);
2647
2648                if (p->Valid(mEpsilon))
2649                {
2650                        candidatePolys.push_back(p);
2651                        candidatePlanes.push_back(*pit);
2652                }
2653        }
2654
2655        // add faces of bounding box (also could be faces of the cell)
2656        for (int i = 0; i < 6; ++ i)
2657        {
2658                VertexContainer vertices;
2659
2660                for (int j = 0; j < 4; ++ j)
2661                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
2662
2663                Polygon3 *poly = new Polygon3(vertices);
2664
2665                candidatePolys.push_back(poly);
2666                candidatePlanes.push_back(poly->GetSupportingPlane());
2667        }
2668
2669        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
2670        {
2671                // polygon is split by all other planes
2672                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
2673                {
2674                        if (i == j) // polygon and plane are coincident
2675                                continue;
2676
2677                        VertexContainer splitPts;
2678                        Polygon3 *frontPoly, *backPoly;
2679
2680                        const int cf =
2681                                candidatePolys[i]->ClassifyPlane(halfSpaces[j],
2682                                                                                                 mEpsilon);
2683
2684                        switch (cf)
2685                        {
2686                                case Polygon3::SPLIT:
2687                                        frontPoly = new Polygon3();
2688                                        backPoly = new Polygon3();
2689
2690                                        candidatePolys[i]->Split(halfSpaces[j],
2691                                                                                         *frontPoly,
2692                                                                                         *backPoly,
2693                                                                                         mEpsilon);
2694
2695                                        DEL_PTR(candidatePolys[i]);
2696
2697                                        if (backPoly->Valid(mEpsilon))
2698                                                candidatePolys[i] = backPoly;
2699                                        else
2700                                                DEL_PTR(backPoly);
2701
2702                                        // outside, don't need this
2703                                        DEL_PTR(frontPoly);
2704                                        break;
2705                                // polygon outside of halfspace
2706                                case Polygon3::FRONT_SIDE:
2707                                        DEL_PTR(candidatePolys[i]);
2708                                        break;
2709                                // just take polygon as it is
2710                                case Polygon3::BACK_SIDE:
2711                                case Polygon3::COINCIDENT:
2712                                default:
2713                                        break;
2714                        }
2715                }
2716
2717                if (candidatePolys[i])
2718                {
2719                        geom.Add(candidatePolys[i], candidatePlanes[i]);
2720                        //      geom.mPolys.push_back(candidates[i]);
2721                }
2722        }
2723}
2724
2725
2726void VspBspTree::ConstructGeometry(ViewCell *vc,
2727                                                                   BspNodeGeometry &vcGeom) const
2728{
2729        ViewCellContainer leaves;
2730       
2731        mViewCellsTree->CollectLeaves(vc, leaves);
2732
2733        ViewCellContainer::const_iterator it, it_end = leaves.end();
2734
2735        for (it = leaves.begin(); it != it_end; ++ it)
2736        {
2737                BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
2738               
2739                ConstructGeometry(l, vcGeom);
2740        }
2741}
2742
2743
2744int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
2745                                                          const bool onlyUnmailed) const
2746{
2747        stack<bspNodePair> nodeStack;
2748       
2749        BspNodeGeometry nodeGeom;
2750        ConstructGeometry(n, nodeGeom);
2751       
2752        // split planes from the root to this node
2753        // needed to verify that we found neighbor leaf
2754        // TODO: really needed?
2755        vector<Plane3> halfSpaces;
2756        ExtractHalfSpaces(n, halfSpaces);
2757
2758
2759        BspNodeGeometry *rgeom = new BspNodeGeometry();
2760        ConstructGeometry(mRoot, *rgeom);
2761
2762        nodeStack.push(bspNodePair(mRoot, rgeom));
2763
2764        while (!nodeStack.empty())
2765        {
2766                BspNode *node = nodeStack.top().first;
2767                BspNodeGeometry *geom = nodeStack.top().second;
2768       
2769                nodeStack.pop();
2770
2771                if (node->IsLeaf())
2772                {
2773                        // test if this leaf is in valid view space
2774                        if (node->TreeValid() &&
2775                                (node != n) &&
2776                                (!onlyUnmailed || !node->Mailed()))
2777                        {
2778                                bool isAdjacent = true;
2779
2780                                if (1)
2781                                {
2782                                        // test all planes of current node if still adjacent
2783                                        for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
2784                                        {
2785                                                const int cf =
2786                                                        Polygon3::ClassifyPlane(geom->GetPolys(),
2787                                                                                                        halfSpaces[i],
2788                                                                                                        mEpsilon);
2789
2790                                                if (cf == Polygon3::FRONT_SIDE)
2791                                                {
2792                                                        isAdjacent = false;
2793                                                }
2794                                        }
2795                                }
2796                                else
2797                                {
2798                                        // TODO: why is this wrong??
2799                                        // test all planes of current node if still adjacent
2800                                        for (int i = 0; (i < nodeGeom.Size()) && isAdjacent; ++ i)
2801                                        {
2802                                                Polygon3 *poly = nodeGeom.GetPolys()[i];
2803
2804                                                const int cf =
2805                                                        Polygon3::ClassifyPlane(geom->GetPolys(),
2806                                                                                                        poly->GetSupportingPlane(),
2807                                                                                                        mEpsilon);
2808
2809                                                if (cf == Polygon3::FRONT_SIDE)
2810                                                {
2811                                                        isAdjacent = false;
2812                                                }
2813                                        }
2814                                }
2815                                // neighbor was found
2816                                if (isAdjacent)
2817                                {       
2818                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2819                                }
2820                        }
2821                }
2822                else
2823                {
2824                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2825
2826                        const int cf = Polygon3::ClassifyPlane(nodeGeom.GetPolys(),
2827                                                                                                   interior->GetPlane(),
2828                                                                                                   mEpsilon);
2829                       
2830                        BspNode *front = interior->GetFront();
2831                        BspNode *back = interior->GetBack();
2832           
2833                        BspNodeGeometry *fGeom = new BspNodeGeometry();
2834                        BspNodeGeometry *bGeom = new BspNodeGeometry();
2835
2836                        geom->SplitGeometry(*fGeom,
2837                                                                *bGeom,
2838                                                                interior->GetPlane(),
2839                                                                mBox,
2840                                                                //0.0000001f);
2841                                                                mEpsilon);
2842               
2843                        if (cf == Polygon3::BACK_SIDE)
2844                        {
2845                                nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2846                                DEL_PTR(fGeom);
2847                        }
2848                        else
2849                        {
2850                                if (cf == Polygon3::FRONT_SIDE)
2851                                {
2852                                        nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2853                                        DEL_PTR(bGeom);
2854                                }
2855                                else
2856                                {       // random decision
2857                                        nodeStack.push(bspNodePair(front, fGeom));
2858                                        nodeStack.push(bspNodePair(back, bGeom));
2859                                }
2860                        }
2861                }
2862       
2863                DEL_PTR(geom);
2864        }
2865
2866        return (int)neighbors.size();
2867}
2868
2869
2870
2871int VspBspTree::FindApproximateNeighbors(BspNode *n,
2872                                                                                 vector<BspLeaf *> &neighbors,
2873                                                                                 const bool onlyUnmailed) const
2874{
2875        stack<bspNodePair> nodeStack;
2876       
2877        BspNodeGeometry nodeGeom;
2878        ConstructGeometry(n, nodeGeom);
2879       
2880        // split planes from the root to this node
2881        // needed to verify that we found neighbor leaf
2882        // TODO: really needed?
2883        vector<Plane3> halfSpaces;
2884        ExtractHalfSpaces(n, halfSpaces);
2885
2886
2887        BspNodeGeometry *rgeom = new BspNodeGeometry();
2888        ConstructGeometry(mRoot, *rgeom);
2889
2890        nodeStack.push(bspNodePair(mRoot, rgeom));
2891
2892        while (!nodeStack.empty())
2893        {
2894                BspNode *node = nodeStack.top().first;
2895                BspNodeGeometry *geom = nodeStack.top().second;
2896       
2897                nodeStack.pop();
2898
2899                if (node->IsLeaf())
2900                {
2901                        // test if this leaf is in valid view space
2902                        if (node->TreeValid() &&
2903                                (node != n) &&
2904                                (!onlyUnmailed || !node->Mailed()))
2905                        {
2906                                bool isAdjacent = true;
2907
2908                                // neighbor was found
2909                                if (isAdjacent)
2910                                {       
2911                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2912                                }
2913                        }
2914                }
2915                else
2916                {
2917                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2918
2919                        const int cf = Polygon3::ClassifyPlane(nodeGeom.GetPolys(),
2920                                                                                                   interior->GetPlane(),
2921                                                                                                   mEpsilon);
2922                       
2923                        BspNode *front = interior->GetFront();
2924                        BspNode *back = interior->GetBack();
2925           
2926                        BspNodeGeometry *fGeom = new BspNodeGeometry();
2927                        BspNodeGeometry *bGeom = new BspNodeGeometry();
2928
2929                        geom->SplitGeometry(*fGeom,
2930                                                                *bGeom,
2931                                                                interior->GetPlane(),
2932                                                                mBox,
2933                                                                //0.0000001f);
2934                                                                mEpsilon);
2935               
2936                        if (cf == Polygon3::BACK_SIDE)
2937                        {
2938                                nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2939                                DEL_PTR(fGeom);
2940                                }
2941                        else
2942                        {
2943                                if (cf == Polygon3::FRONT_SIDE)
2944                                {
2945                                        nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2946                                        DEL_PTR(bGeom);
2947                                }
2948                                else
2949                                {       // random decision
2950                                        nodeStack.push(bspNodePair(front, fGeom));
2951                                        nodeStack.push(bspNodePair(back, bGeom));
2952                                }
2953                        }
2954                }
2955       
2956                DEL_PTR(geom);
2957        }
2958
2959        return (int)neighbors.size();
2960}
2961
2962
2963
2964BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
2965{
2966    stack<BspNode *> nodeStack;
2967        nodeStack.push(mRoot);
2968
2969        int mask = rand();
2970
2971        while (!nodeStack.empty())
2972        {
2973                BspNode *node = nodeStack.top();
2974                nodeStack.pop();
2975
2976                if (node->IsLeaf())
2977                {
2978                        return dynamic_cast<BspLeaf *>(node);
2979                }
2980                else
2981                {
2982                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2983                        BspNode *next;
2984                        BspNodeGeometry geom;
2985
2986                        // todo: not very efficient: constructs full cell everytime
2987                        ConstructGeometry(interior, geom);
2988
2989                        const int cf =
2990                                Polygon3::ClassifyPlane(geom.GetPolys(), halfspace, mEpsilon);
2991
2992                        if (cf == Polygon3::BACK_SIDE)
2993                                next = interior->GetFront();
2994                        else
2995                                if (cf == Polygon3::FRONT_SIDE)
2996                                        next = interior->GetFront();
2997                        else
2998                        {
2999                                // random decision
3000                                if (mask & 1)
3001                                        next = interior->GetBack();
3002                                else
3003                                        next = interior->GetFront();
3004                                mask = mask >> 1;
3005                        }
3006
3007                        nodeStack.push(next);
3008                }
3009        }
3010
3011        return NULL;
3012}
3013
3014
3015BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
3016{
3017        stack<BspNode *> nodeStack;
3018
3019        nodeStack.push(mRoot);
3020
3021        int mask = rand();
3022
3023        while (!nodeStack.empty())
3024        {
3025                BspNode *node = nodeStack.top();
3026                nodeStack.pop();
3027
3028                if (node->IsLeaf())
3029                {
3030                        if ( (!onlyUnmailed || !node->Mailed()) )
3031                                return dynamic_cast<BspLeaf *>(node);
3032                }
3033                else
3034                {
3035                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
3036
3037                        // random decision
3038                        if (mask & 1)
3039                                nodeStack.push(interior->GetBack());
3040                        else
3041                                nodeStack.push(interior->GetFront());
3042
3043                        mask = mask >> 1;
3044                }
3045        }
3046
3047        return NULL;
3048}
3049
3050
3051int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
3052{
3053        int pvsSize = 0;
3054
3055        RayInfoContainer::const_iterator rit, rit_end = rays.end();
3056
3057        Intersectable::NewMail();
3058
3059        for (rit = rays.begin(); rit != rays.end(); ++ rit)
3060        {
3061                VssRay *ray = (*rit).mRay;
3062
3063                if (ray->mOriginObject)
3064                {
3065                        if (!ray->mOriginObject->Mailed())
3066                        {
3067                                ray->mOriginObject->Mail();
3068                                ++ pvsSize;
3069                        }
3070                }
3071                if (ray->mTerminationObject)
3072                {
3073                        if (!ray->mTerminationObject->Mailed())
3074                        {
3075                                ray->mTerminationObject->Mail();
3076                                ++ pvsSize;
3077                        }
3078                }
3079        }
3080
3081        return pvsSize;
3082}
3083
3084
3085float VspBspTree::GetEpsilon() const
3086{
3087        return mEpsilon;
3088}
3089
3090
3091int VspBspTree::SplitPolygons(const Plane3 &plane,
3092                                                          PolygonContainer &polys,
3093                                                          PolygonContainer &frontPolys,
3094                                                          PolygonContainer &backPolys,
3095                                                          PolygonContainer &coincident) const
3096{
3097        int splits = 0;
3098
3099        PolygonContainer::const_iterator it, it_end = polys.end();
3100
3101        for (it = polys.begin(); it != polys.end(); ++ it)     
3102        {
3103                Polygon3 *poly = *it;
3104
3105                // classify polygon
3106                const int cf = poly->ClassifyPlane(plane, mEpsilon);
3107
3108                switch (cf)
3109                {
3110                        case Polygon3::COINCIDENT:
3111                                coincident.push_back(poly);
3112                                break;
3113                        case Polygon3::FRONT_SIDE:
3114                                frontPolys.push_back(poly);
3115                                break;
3116                        case Polygon3::BACK_SIDE:
3117                                backPolys.push_back(poly);
3118                                break;
3119                        case Polygon3::SPLIT:
3120                                backPolys.push_back(poly);
3121                                frontPolys.push_back(poly);
3122                                ++ splits;
3123                                break;
3124                        default:
3125                Debug << "SHOULD NEVER COME HERE\n";
3126                                break;
3127                }
3128        }
3129
3130        return splits;
3131}
3132
3133
3134int VspBspTree::CastLineSegment(const Vector3 &origin,
3135                                                                const Vector3 &termination,
3136                                                                vector<ViewCell *> &viewcells)
3137{
3138        int hits = 0;
3139        stack<BspRayTraversalData> tStack;
3140
3141        float mint = 0.0f, maxt = 1.0f;
3142
3143        Intersectable::NewMail();
3144        ViewCell::NewMail();
3145
3146        Vector3 entp = origin;
3147        Vector3 extp = termination;
3148
3149        BspNode *node = mRoot;
3150        BspNode *farChild = NULL;
3151
3152        float t;
3153        const float thresh = 1e-6f; // matt: change this
3154       
3155        while (1)
3156        {
3157                if (!node->IsLeaf())
3158                {
3159                        BspInterior *in = dynamic_cast<BspInterior *>(node);
3160
3161                        Plane3 splitPlane = in->GetPlane();
3162                       
3163                        const int entSide = splitPlane.Side(entp, thresh);
3164                        const int extSide = splitPlane.Side(extp, thresh);
3165
3166                        if (entSide < 0)
3167                        {
3168                                node = in->GetBack();
3169                               
3170                                // plane does not split ray => no far child
3171                                if (extSide <= 0)
3172                                        continue;
3173 
3174                                farChild = in->GetFront(); // plane splits ray
3175                        }
3176                        else if (entSide > 0)
3177                        {
3178                                node = in->GetFront();
3179
3180                                if (extSide >= 0) // plane does not split ray => no far child
3181                                        continue;
3182
3183                                farChild = in->GetBack(); // plane splits ray
3184                        }
3185                        else // one of the ray end points is on the plane
3186                        {       // NOTE: what to do if ray is coincident with plane?
3187                                if (extSide < 0)
3188                                        node = in->GetBack();
3189                                else //if (extSide > 0)
3190                                        node = in->GetFront();
3191                                //else break; // coincident => count no intersections
3192
3193                                continue; // no far child
3194                        }
3195
3196                        // push data for far child
3197                        tStack.push(BspRayTraversalData(farChild, extp));
3198
3199                        // find intersection of ray segment with plane
3200                        extp = splitPlane.FindIntersection(origin, extp, &t);
3201                }
3202                else
3203                {
3204                        // reached leaf => intersection with view cell
3205                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
3206                        ViewCell *viewCell;
3207                       
3208                        if (0)
3209                                viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
3210                        else
3211                                viewCell = leaf->GetViewCell();
3212
3213                        if (!viewCell->Mailed())
3214                        {
3215                                viewcells.push_back(viewCell);
3216                                viewCell->Mail();
3217                                ++ hits;
3218                        }
3219
3220                        //-- fetch the next far child from the stack
3221                        if (tStack.empty())
3222                                break;
3223
3224                        entp = extp;
3225                       
3226                        const BspRayTraversalData &s = tStack.top();
3227
3228                        node = s.mNode;
3229                        extp = s.mExitPoint;
3230
3231                        tStack.pop();
3232                }
3233        }
3234
3235        return hits;
3236}
3237
3238
3239
3240
3241int VspBspTree::TreeDistance(BspNode *n1, BspNode *n2) const
3242{
3243        std::deque<BspNode *> path1;
3244        BspNode *p1 = n1;
3245
3246        // create path from node 1 to root
3247        while (p1)
3248        {
3249                if (p1 == n2) // second node on path
3250                        return (int)path1.size();
3251
3252                path1.push_front(p1);
3253                p1 = p1->GetParent();
3254        }
3255
3256        int depth = n2->GetDepth();
3257        int d = depth;
3258
3259        BspNode *p2 = n2;
3260
3261        // compare with same depth
3262        while (1)
3263        {
3264                if ((d < (int)path1.size()) && (p2 == path1[d]))
3265                        return (depth - d) + ((int)path1.size() - 1 - d);
3266
3267                -- d;
3268                p2 = p2->GetParent();
3269        }
3270
3271        return 0; // never come here
3272}
3273
3274
3275BspNode *VspBspTree::CollapseTree(BspNode *node, int &collapsed)
3276{
3277// TODO
3278#if HAS_TO_BE_REDONE
3279        if (node->IsLeaf())
3280                return node;
3281
3282        BspInterior *interior = dynamic_cast<BspInterior *>(node);
3283
3284        BspNode *front = CollapseTree(interior->GetFront(), collapsed);
3285        BspNode *back = CollapseTree(interior->GetBack(), collapsed);
3286
3287        if (front->IsLeaf() && back->IsLeaf())
3288        {
3289                BspLeaf *frontLeaf = dynamic_cast<BspLeaf *>(front);
3290                BspLeaf *backLeaf = dynamic_cast<BspLeaf *>(back);
3291
3292                //-- collapse tree
3293                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
3294                {
3295                        BspViewCell *vc = frontLeaf->GetViewCell();
3296
3297                        BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);
3298                        leaf->SetTreeValid(frontLeaf->TreeValid());
3299
3300                        // replace a link from node's parent
3301                        if (leaf->GetParent())
3302                                leaf->GetParent()->ReplaceChildLink(node, leaf);
3303                        else
3304                                mRoot = leaf;
3305
3306                        ++ collapsed;
3307                        delete interior;
3308
3309                        return leaf;
3310                }
3311        }
3312#endif
3313        return node;
3314}
3315
3316
3317int VspBspTree::CollapseTree()
3318{
3319        int collapsed = 0;
3320        //TODO
3321#if HAS_TO_BE_REDONE
3322        (void)CollapseTree(mRoot, collapsed);
3323
3324        // revalidate leaves
3325        RepairViewCellsLeafLists();
3326#endif
3327        return collapsed;
3328}
3329
3330
3331void VspBspTree::RepairViewCellsLeafLists()
3332{
3333// TODO
3334#if HAS_TO_BE_REDONE
3335        // list not valid anymore => clear
3336        stack<BspNode *> nodeStack;
3337        nodeStack.push(mRoot);
3338
3339        ViewCell::NewMail();
3340
3341        while (!nodeStack.empty())
3342        {
3343                BspNode *node = nodeStack.top();
3344                nodeStack.pop();
3345
3346                if (node->IsLeaf())
3347                {
3348                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
3349
3350                        BspViewCell *viewCell = leaf->GetViewCell();
3351
3352                        if (!viewCell->Mailed())
3353                        {
3354                                viewCell->mLeaves.clear();
3355                                viewCell->Mail();
3356                        }
3357       
3358                        viewCell->mLeaves.push_back(leaf);
3359
3360                }
3361                else
3362                {
3363                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
3364
3365                        nodeStack.push(interior->GetFront());
3366                        nodeStack.push(interior->GetBack());
3367                }
3368        }
3369// TODO
3370#endif
3371}
3372
3373
3374typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;
3375
3376
3377int VspBspTree::CastBeam(Beam &beam)
3378{
3379    stack<bspNodePair> nodeStack;
3380        BspNodeGeometry *rgeom = new BspNodeGeometry();
3381        ConstructGeometry(mRoot, *rgeom);
3382
3383        nodeStack.push(bspNodePair(mRoot, rgeom));
3384 
3385        ViewCell::NewMail();
3386
3387        while (!nodeStack.empty())
3388        {
3389                BspNode *node = nodeStack.top().first;
3390                BspNodeGeometry *geom = nodeStack.top().second;
3391                nodeStack.pop();
3392               
3393                AxisAlignedBox3 box;
3394                geom->GetBoundingBox(box);
3395
3396                const int side = beam.ComputeIntersection(box);
3397               
3398                switch (side)
3399                {
3400                case -1:
3401                        CollectViewCells(node, true, beam.mViewCells, true);
3402                        break;
3403                case 0:
3404                       
3405                        if (node->IsLeaf())
3406                        {
3407                                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
3408                       
3409                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
3410                                {
3411                                        leaf->GetViewCell()->Mail();
3412                                        beam.mViewCells.push_back(leaf->GetViewCell());
3413                                }
3414                        }
3415                        else
3416                        {
3417                                BspInterior *interior = dynamic_cast<BspInterior *>(node);
3418                       
3419                                BspNode *first = interior->GetFront();
3420                                BspNode *second = interior->GetBack();
3421           
3422                                BspNodeGeometry *firstGeom = new BspNodeGeometry();
3423                                BspNodeGeometry *secondGeom = new BspNodeGeometry();
3424
3425                                geom->SplitGeometry(*firstGeom,
3426                                                                        *secondGeom,
3427                                                                        interior->GetPlane(),
3428                                                                        mBox,
3429                                                                        //0.0000001f);
3430                                                                        mEpsilon);
3431
3432                                // decide on the order of the nodes
3433                                if (DotProd(beam.mPlanes[0].mNormal,
3434                                        interior->GetPlane().mNormal) > 0)
3435                                {
3436                                        swap(first, second);
3437                                        swap(firstGeom, secondGeom);
3438                                }
3439
3440                                nodeStack.push(bspNodePair(first, firstGeom));
3441                                nodeStack.push(bspNodePair(second, secondGeom));
3442                        }
3443                       
3444                        break;
3445                default:
3446                        // default: cull
3447                        break;
3448                }
3449               
3450                DEL_PTR(geom);
3451               
3452        }
3453
3454        return (int)beam.mViewCells.size();
3455}
3456
3457
3458void VspBspTree::SetViewCellsManager(ViewCellsManager *vcm)
3459{
3460        mViewCellsManager = vcm;
3461}
3462
3463
3464int VspBspTree::CollectMergeCandidates(const vector<BspLeaf *> leaves,
3465                                                                           vector<MergeCandidate> &candidates)
3466{
3467        BspLeaf::NewMail();
3468       
3469        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
3470
3471        int numCandidates = 0;
3472
3473        // find merge candidates and push them into queue
3474        for (it = leaves.begin(); it != it_end; ++ it)
3475        {
3476                BspLeaf *leaf = *it;
3477               
3478                // the same leaves must not be part of two merge candidates
3479                leaf->Mail();
3480               
3481                vector<BspLeaf *> neighbors;
3482               
3483                // appoximate neighbor search has slightl relaxed constraints
3484                if (1)
3485                        FindNeighbors(leaf, neighbors, true);
3486                else
3487                        FindApproximateNeighbors(leaf, neighbors, true);
3488
3489                vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();
3490
3491                // TODO: test if at least one ray goes from one leaf to the other
3492                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
3493                {
3494                        if ((*nit)->GetViewCell() != leaf->GetViewCell())
3495                        {
3496                                MergeCandidate mc(leaf->GetViewCell(), (*nit)->GetViewCell());
3497
3498                                // dont't merge view cells if they are empty and not front and back leaf of the same parent
3499                                // in the tree?
3500                                if (mEmptyViewCellsMergeAllowed ||
3501                                        !leaf->GetViewCell()->GetPvs().Empty() ||
3502                                        !(*nit)->GetViewCell()->GetPvs().Empty() ||
3503                    leaf->IsSibling(*nit))
3504                                {
3505                                        candidates.push_back(mc);
3506                                }
3507
3508                                ++ numCandidates;
3509                                if ((numCandidates % 1000) == 0)
3510                                {
3511                                        cout << "collected " << numCandidates << " merge candidates" << endl;
3512                                }
3513                        }
3514                }
3515        }
3516
3517        Debug << "merge queue: " << (int)candidates.size() << endl;
3518        Debug << "leaves in queue: " << numCandidates << endl;
3519       
3520
3521        return (int)leaves.size();
3522}
3523
3524
3525int VspBspTree::CollectMergeCandidates(const VssRayContainer &rays,
3526                                                                           vector<MergeCandidate> &candidates)
3527{
3528        ViewCell::NewMail();
3529        long startTime = GetTime();
3530       
3531        map<BspLeaf *, vector<BspLeaf*> > neighborMap;
3532        ViewCellContainer::const_iterator iit;
3533
3534        int numLeaves = 0;
3535       
3536        BspLeaf::NewMail();
3537
3538        for (int i = 0; i < (int)rays.size(); ++ i)
3539        { 
3540                VssRay *ray = rays[i];
3541       
3542                // traverse leaves stored in the rays and compare and
3543                // merge consecutive leaves (i.e., the neighbors in the tree)
3544                if (ray->mViewCells.size() < 2)
3545                        continue;
3546//TODO viewcellhierarchy
3547                iit = ray->mViewCells.begin();
3548                BspViewCell *bspVc = dynamic_cast<BspViewCell *>(*(iit ++));
3549                BspLeaf *leaf = bspVc->mLeaf;
3550               
3551                // traverse intersections
3552                // consecutive leaves are neighbors => add them to queue
3553                for (; iit != ray->mViewCells.end(); ++ iit)
3554                {
3555                        // next pair
3556                        BspLeaf *prevLeaf = leaf;
3557                        bspVc = dynamic_cast<BspViewCell *>(*iit);
3558            leaf = bspVc->mLeaf;
3559
3560                        // view space not valid or same view cell
3561                        if (!leaf->TreeValid() || !prevLeaf->TreeValid() ||
3562                                (leaf->GetViewCell() == prevLeaf->GetViewCell()))
3563                                continue;
3564
3565                vector<BspLeaf *> &neighbors = neighborMap[leaf];
3566                       
3567                        bool found = false;
3568
3569                        // both leaves inserted in queue already =>
3570                        // look if double pair already exists
3571                        if (leaf->Mailed() && prevLeaf->Mailed())
3572                        {
3573                                vector<BspLeaf *>::const_iterator it, it_end = neighbors.end();
3574                               
3575                for (it = neighbors.begin(); !found && (it != it_end); ++ it)
3576                                        if (*it == prevLeaf)
3577                                                found = true; // already in queue
3578                        }
3579               
3580                        if (!found)
3581                        {
3582                                // this pair is not in map yet
3583                                // => insert into the neighbor map and the queue
3584                                neighbors.push_back(prevLeaf);
3585                                neighborMap[prevLeaf].push_back(leaf);
3586
3587                                leaf->Mail();
3588                                prevLeaf->Mail();
3589               
3590                                MergeCandidate mc(leaf->GetViewCell(), prevLeaf->GetViewCell());
3591                               
3592                                candidates.push_back(mc);
3593
3594                                if (((int)candidates.size() % 1000) == 0)
3595                                {
3596                                        cout << "collected " << (int)candidates.size() << " merge candidates" << endl;
3597                                }
3598                        }
3599        }
3600        }
3601
3602        Debug << "neighbormap size: " << (int)neighborMap.size() << endl;
3603        Debug << "merge queue: " << (int)candidates.size() << endl;
3604        Debug << "leaves in queue: " << numLeaves << endl;
3605
3606
3607        //-- collect the leaves which haven't been found by ray casting
3608        if (0)
3609        {
3610                cout << "finding additional merge candidates using geometry" << endl;
3611                vector<BspLeaf *> leaves;
3612                CollectLeaves(leaves, true);
3613                Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
3614                CollectMergeCandidates(leaves, candidates);
3615        }
3616
3617        return numLeaves;
3618}
3619
3620
3621
3622
3623ViewCell *VspBspTree::GetViewCell(const Vector3 &point)
3624{
3625  if (mRoot == NULL)
3626        return NULL;
3627 
3628  stack<BspNode *> nodeStack;
3629  nodeStack.push(mRoot);
3630 
3631  ViewCell *viewcell = NULL;
3632 
3633  while (!nodeStack.empty())  {
3634        BspNode *node = nodeStack.top();
3635        nodeStack.pop();
3636       
3637        if (node->IsLeaf()) {
3638          viewcell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
3639          break;
3640        } else {
3641         
3642          BspInterior *interior = dynamic_cast<BspInterior *>(node);
3643               
3644          // random decision
3645          if (interior->GetPlane().Side(point) < 0)
3646                nodeStack.push(interior->GetBack());
3647          else
3648                nodeStack.push(interior->GetFront());
3649        }
3650  }
3651 
3652  return viewcell;
3653}
3654
3655
3656bool VspBspTree::ViewPointValid(const Vector3 &viewPoint) const
3657{
3658        BspNode *node = mRoot;
3659
3660        while (1)
3661        {
3662                // early exit
3663                if (node->TreeValid())
3664                        return true;
3665
3666                if (node->IsLeaf())
3667                        return false;
3668                       
3669                BspInterior *in = dynamic_cast<BspInterior *>(node);
3670                                       
3671                if (in->GetPlane().Side(viewPoint) <= 0)
3672                {
3673                        node = in->GetBack();
3674                }
3675                else
3676                {
3677                        node = in->GetFront();
3678                }
3679        }
3680
3681        // should never come here
3682        return false;
3683}
3684
3685
3686void VspBspTree::PropagateUpValidity(BspNode *node)
3687{
3688        const bool isValid = node->TreeValid();
3689
3690        // propagative up invalid flag until only invalid nodes exist over this node
3691        if (!isValid)
3692        {
3693                while (!node->IsRoot() && node->GetParent()->TreeValid())
3694                {
3695                        node = node->GetParent();
3696                        node->SetTreeValid(false);
3697                }
3698        }
3699        else
3700        {
3701                // propagative up valid flag until one of the subtrees is invalid
3702                while (!node->IsRoot() && !node->TreeValid())
3703                {
3704            node = node->GetParent();
3705                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
3706                       
3707                        // the parent is valid iff both leaves are valid
3708                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
3709                                                           interior->GetFront()->TreeValid());
3710                }
3711        }
3712}
3713
3714
3715bool VspBspTree::Export(ofstream &stream)
3716{
3717        ExportNode(mRoot, stream);
3718
3719        return true;
3720}
3721
3722
3723void VspBspTree::ExportNode(BspNode *node, ofstream &stream)
3724{
3725        if (node->IsLeaf())
3726        {
3727                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
3728                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
3729
3730                int id = -1;
3731                if (viewCell != mOutOfBoundsCell)
3732                        id = viewCell->GetId();
3733
3734                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
3735        }
3736        else
3737        {
3738                BspInterior *interior = dynamic_cast<BspInterior *>(node);
3739       
3740                Plane3 plane = interior->GetPlane();
3741                stream << "<Interior plane=\"" << plane.mNormal.x << " "
3742                           << plane.mNormal.y << " " << plane.mNormal.z << " "
3743                           << plane.mD << "\">" << endl;
3744
3745                ExportNode(interior->GetBack(), stream);
3746                ExportNode(interior->GetFront(), stream);
3747
3748                stream << "</Interior>" << endl;
3749        }
3750}
Note: See TracBrowser for help on using the repository browser.