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

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