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

Revision 810, 92.2 KB checked in by mattausch, 18 years ago (diff)

added from point visibility tree:
this tree is constructed using from point visibility information

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