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

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