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

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