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

Revision 648, 75.9 KB checked in by mattausch, 18 years ago (diff)

removed problems with bsp due to polygons equal to scene box faces. removing
this polyogns in a preprocessing step. also using this for vsp bsp tree, because it
is generally a good thing.

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