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

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