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

Revision 650, 76.0 KB checked in by mattausch, 19 years ago (diff)

added spatial hierarchy based merge tree

Line 
1#include <stack>
2#include <time.h>
3#include <iomanip>
4
5#include "Plane3.h"
6#include "VspBspTree.h"
7#include "Mesh.h"
8#include "common.h"
9#include "ViewCell.h"
10#include "Environment.h"
11#include "Polygon3.h"
12#include "Ray.h"
13#include "AxisAlignedBox3.h"
14#include "Exporter.h"
15#include "Plane3.h"
16#include "ViewCellBsp.h"
17#include "ViewCellsManager.h"
18#include "Beam.h"
19
20#define USE_FIXEDPOINT_T 0
21
22
23//-- static members
24
25/** Evaluates split plane classification with respect to the plane's
26        contribution for a minimum number of ray splits.
27*/
28const float VspBspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
29/** Evaluates split plane classification with respect to the plane's
30        contribution for  balanced rays.
31*/
32const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
33
34
35int VspBspTree::sFrontId = 0;
36int VspBspTree::sBackId = 0;
37int VspBspTree::sFrontAndBackId = 0;
38
39
40
41
42// pvs penalty can be different from pvs size
43inline float EvalPvsPenalty(const int pvs,
44                                                        const int lower,
45                                                        const int upper)
46{
47        // clamp to minmax values
48        if (pvs < lower)
49                return (float)lower;
50        if (pvs > upper)
51                return (float)upper;
52
53        return (float)pvs;
54}
55
56
57
58
59/******************************************************************************/
60/*                       class VspBspTree implementation                      */
61/******************************************************************************/
62
63
64VspBspTree::VspBspTree():
65mRoot(NULL),
66mUseAreaForPvs(false),
67mCostNormalizer(Limits::Small),
68mViewCellsManager(NULL),
69mOutOfBoundsCell(NULL),
70mStoreRays(false),
71mRenderCostWeight(0.5),
72mUseRandomAxis(false),
73mTimeStamp(1)
74{
75        bool randomize = false;
76        environment->GetBoolValue("VspBspTree.Construction.randomize", randomize);
77        if (randomize)
78                Randomize(); // initialise random generator for heuristics
79
80        //-- termination criteria for autopartition
81        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
82        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
83        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
84        environment->GetFloatValue("VspBspTree.Termination.minProbability", mTermMinProbability);
85        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
86        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
87        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
88        environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
89        environment->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
90        //-- max cost ratio for early tree termination
91        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
92
93        // 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" << (float)mTotalPvsSize / (float)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        interior->mTimeStamp = mTimeStamp ++;
769        //frontData.mNode->mTimeStamp = ++ mTimeStamp;
770        //backData.mNode->mTimeStamp = mTimeStamp;
771
772        //DEL_PTR(leaf);
773        return interior;
774}
775
776
777void VspBspTree::AddToPvs(BspLeaf *leaf,
778                                                  const RayInfoContainer &rays,
779                                                  float &sampleContributions,
780                                                  int &contributingSamples)
781{
782  sampleContributions = 0;
783  contributingSamples = 0;
784 
785  RayInfoContainer::const_iterator it, it_end = rays.end();
786 
787  ViewCell *vc = leaf->GetViewCell();
788 
789  // add contributions from samples to the PVS
790  for (it = rays.begin(); it != it_end; ++ it)
791        {
792          float sc = 0.0f;
793          VssRay *ray = (*it).mRay;
794          bool madeContrib = false;
795          float contribution;
796          if (ray->mTerminationObject) {
797                if (vc->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf, contribution))
798                  madeContrib = true;
799                sc += contribution;
800          }
801         
802          if (ray->mOriginObject) {
803                if (vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf, contribution))
804                  madeContrib = true;
805                sc += contribution;
806          }
807         
808          sampleContributions += sc;
809          if (madeContrib)
810                  ++ contributingSamples;
811               
812          //leaf->mVssRays.push_back(ray);
813        }
814}
815
816
817void VspBspTree::SortSplitCandidates(const RayInfoContainer &rays, const int axis)
818{
819        mSplitCandidates->clear();
820
821        int requestedSize = 2 * (int)(rays.size());
822        // creates a sorted split candidates array
823        if (mSplitCandidates->capacity() > 500000 &&
824                requestedSize < (int)(mSplitCandidates->capacity()  / 10) )
825        {
826        delete mSplitCandidates;
827                mSplitCandidates = new vector<SortableEntry>;
828        }
829
830        mSplitCandidates->reserve(requestedSize);
831
832        // insert all queries
833        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
834        {
835                const bool positive = (*ri).mRay->HasPosDir(axis);
836                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
837                                                                                                  (*ri).ExtrapOrigin(axis), (*ri).mRay));
838
839                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
840                                                                                                  (*ri).ExtrapTermination(axis), (*ri).mRay));
841        }
842
843        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
844}
845
846
847float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
848                                                                                  const AxisAlignedBox3 &box,
849                                                                                  const int pvsSize,
850                                                                                  const int &axis,
851                                          float &position)
852{
853        SortSplitCandidates(rays, axis);
854
855        // go through the lists, count the number of objects left and right
856        // and evaluate the following cost funcion:
857        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
858
859        int pvsl = 0, pvsr = pvsSize;
860
861        int pvsBack = pvsl;
862        int pvsFront = pvsr;
863
864        float minBox = box.Min(axis);
865        float maxBox = box.Max(axis);
866        float sizeBox = maxBox - minBox;
867
868        float minBand = minBox + 0.1f * (maxBox - minBox);
869        float maxBand = minBox + 0.9f * (maxBox - minBox);
870
871        float sum = (float)pvsSize * sizeBox;
872        float minSum = 1e20f;
873
874        Intersectable::NewMail();
875
876        RayInfoContainer::const_iterator ri, ri_end = rays.end();
877
878        // set all object as belonging to the front pvs
879        for(ri = rays.begin(); ri != ri_end; ++ ri)
880        {
881                Intersectable *oObject = (*ri).mRay->mOriginObject;
882                Intersectable *tObject = (*ri).mRay->mTerminationObject;
883
884                if (oObject)
885                {
886                        if (!oObject->Mailed())
887                        {
888                                oObject->Mail();
889                                oObject->mCounter = 1;
890                        }
891                        else
892                        {
893                                ++ oObject->mCounter;
894                        }
895                }
896                if (tObject)
897                {
898                        if (!tObject->Mailed())
899                        {
900                                tObject->Mail();
901                                tObject->mCounter = 1;
902                        }
903                        else
904                        {
905                                ++ tObject->mCounter;
906                        }
907                }
908        }
909
910        Intersectable::NewMail();
911
912        vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
913
914        for (ci = mSplitCandidates->begin(); ci < ci_end; ++ ci)
915        {
916                VssRay *ray;
917                ray = (*ci).ray;
918               
919                Intersectable *oObject = ray->mOriginObject;
920                Intersectable *tObject = ray->mTerminationObject;
921               
922
923                switch ((*ci).type)
924                {
925                        case SortableEntry::ERayMin:
926                                {
927                                        if (oObject && !oObject->Mailed())
928                                        {
929                                                oObject->Mail();
930                                                ++ pvsl;
931                                        }
932                                        if (tObject && !tObject->Mailed())
933                                        {
934                                                tObject->Mail();
935                                                ++ pvsl;
936                                        }
937                                        break;
938                                }
939                        case SortableEntry::ERayMax:
940                                {
941                                        if (oObject)
942                                        {
943                                                if (-- oObject->mCounter == 0)
944                                                        -- pvsr;
945                                        }
946
947                                        if (tObject)
948                                        {
949                                                if (-- tObject->mCounter == 0)
950                                                        -- pvsr;
951                                        }
952
953                                        break;
954                                }
955                }
956
957                // Note: sufficient to compare size of bounding boxes of front and back side?
958                if ((*ci).value > minBand && (*ci).value < maxBand)
959                {
960                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
961
962                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
963                        // cout<<"cost= "<<sum<<endl;
964
965                        if (sum < minSum)
966                        {
967                                minSum = sum;
968                                position = (*ci).value;
969
970                                pvsBack = pvsl;
971                                pvsFront = pvsr;
972                        }
973                }
974        }
975
976        // -- compute cost
977
978        const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
979        const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
980
981        const float pOverall = sizeBox;
982
983        const float pBack = position - minBox;
984        const float pFront = maxBox - position;
985       
986        const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
987    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
988        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
989       
990        const float oldRenderCost = penaltyOld * pOverall;
991        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
992
993        float ratio = mPvsFactor * newRenderCost / (oldRenderCost + Limits::Small);
994
995        //Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
996        //     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
997
998        return ratio;
999}
1000
1001
1002float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
1003                                                                                 const VspBspTraversalData &tData,
1004                                                                                 int &axis,
1005                                                                                 BspNodeGeometry **frontGeom,
1006                                                                                 BspNodeGeometry **backGeom,
1007                                                                                 float &pFront,
1008                                                                                 float &pBack,
1009                                                                                 const bool useKdSplit)
1010{
1011        float nPosition[3];
1012        float nCostRatio[3];
1013        float nProbFront[3];
1014        float nProbBack[3];
1015
1016        BspNodeGeometry *nFrontGeom[3];
1017        BspNodeGeometry *nBackGeom[3];
1018
1019        for (int i = 0; i < 3; ++ i)
1020        {
1021                nFrontGeom[i] = NULL;
1022                nBackGeom[i] = NULL;
1023        }
1024
1025        int bestAxis = -1;
1026
1027        // create bounding box of node geometry
1028        AxisAlignedBox3 box;
1029        box.Initialize();
1030       
1031        //TODO: for kd split geometry already is box => only take minmax vertices
1032        if (1)
1033        {
1034                PolygonContainer::const_iterator it, it_end = tData.mGeometry->mPolys.end();
1035
1036                for(it = tData.mGeometry->mPolys.begin(); it < it_end; ++ it)
1037                        box.Include(*(*it));
1038        }
1039        else
1040        {
1041                RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();
1042
1043                for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
1044                        box.Include((*ri).ExtrapTermination());
1045        }
1046#if OCTREE_HACK
1047        //Debug << "choosing axis:" << tData.mAxis << endl;
1048        const int sAxis = tData.mAxis;
1049#else
1050        const int sAxis = mUseRandomAxis ? Random(3) : box.Size().DrivingAxis();
1051#endif
1052        for (axis = 0; axis < 3; ++ axis)
1053        {
1054                if (!mOnlyDrivingAxis || (axis == sAxis))
1055                {
1056                        if (!mUseCostHeuristics)
1057                        {
1058                                nFrontGeom[axis] = new BspNodeGeometry();
1059                                nBackGeom[axis] = new BspNodeGeometry();
1060
1061                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1062                                Vector3 normal(0,0,0); normal[axis] = 1.0f;
1063
1064                                // allows faster split because we have axis aligned kd tree boxes
1065                                if (useKdSplit)
1066                                {
1067                                        nCostRatio[axis] = EvalAxisAlignedSplitCost(tData,
1068                                                                                                                                box,
1069                                                                                                                                axis,
1070                                                                                                                                nPosition[axis],
1071                                                                                                                                nProbFront[axis],
1072                                                                                                                                nProbBack[axis]);
1073                                       
1074                                        Vector3 pos;
1075                                       
1076                                        pos = box.Max(); pos[axis] = nPosition[axis];
1077                                        AxisAlignedBox3 bBox(box.Min(), pos);
1078                                        bBox.ExtractPolys(nBackGeom[axis]->mPolys);
1079                                       
1080                                        pos = box.Min(); pos[axis] = nPosition[axis];
1081                                        AxisAlignedBox3 fBox(pos, box.Max());
1082                                        fBox.ExtractPolys(nFrontGeom[axis]->mPolys);
1083                                }
1084                                else
1085                                {
1086                                        nCostRatio[axis] =
1087                                                EvalSplitPlaneCost(Plane3(normal, nPosition[axis]),
1088                                                                                   tData, *nFrontGeom[axis], *nBackGeom[axis],
1089                                                                                   nProbFront[axis], nProbBack[axis]);
1090                                }
1091                        }
1092                        else
1093                        {
1094                                nCostRatio[axis] =
1095                                        BestCostRatioHeuristics(*tData.mRays,
1096                                                                                    box,
1097                                                                                        tData.mPvs,
1098                                                                                        axis,
1099                                                                                        nPosition[axis]);
1100                        }
1101
1102                        if (bestAxis == -1)
1103                        {
1104                                bestAxis = axis;
1105                        }
1106                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
1107                        {
1108                                bestAxis = axis;
1109                        }
1110
1111                }
1112        }
1113
1114        //-- assign values
1115        axis = bestAxis;
1116        pFront = nProbFront[bestAxis];
1117        pBack = nProbBack[bestAxis];
1118
1119        // assign best split nodes geometry
1120        *frontGeom = nFrontGeom[bestAxis];
1121        *backGeom = nBackGeom[bestAxis];
1122
1123        // and delete other geometry
1124        DEL_PTR(nFrontGeom[(bestAxis + 1) % 3]);
1125        DEL_PTR(nBackGeom[(bestAxis + 2) % 3]);
1126
1127        //-- split plane
1128    Vector3 normal(0,0,0); normal[bestAxis] = 1;
1129        plane = Plane3(normal, nPosition[bestAxis]);
1130
1131        return nCostRatio[bestAxis];
1132}
1133
1134
1135bool VspBspTree::SelectPlane(Plane3 &bestPlane,
1136                                                         BspLeaf *leaf,
1137                                                         VspBspTraversalData &data,
1138                                                         VspBspTraversalData &frontData,
1139                                                         VspBspTraversalData &backData,
1140                                                         int &splitAxis)
1141{
1142        // simplest strategy: just take next polygon
1143        if (mSplitPlaneStrategy & RANDOM_POLYGON)
1144        {
1145        if (!data.mPolygons->empty())
1146                {
1147                        const int randIdx =
1148                                (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
1149                        Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
1150
1151                        bestPlane = nextPoly->GetSupportingPlane();
1152                        return true;
1153                }
1154        }
1155
1156        //-- use heuristics to find appropriate plane
1157
1158        // intermediate plane
1159        Plane3 plane;
1160        float lowestCost = MAX_FLOAT;
1161       
1162        // decides if the first few splits should be only axisAligned
1163        const bool onlyAxisAligned  =
1164                (mSplitPlaneStrategy & AXIS_ALIGNED) &&
1165                ((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
1166                ((int)data.GetAvgRayContribution() < mTermMaxRayContriForAxisAligned);
1167       
1168        const int limit = onlyAxisAligned ? 0 :
1169                Min((int)data.mPolygons->size(), mMaxPolyCandidates);
1170
1171        float candidateCost;
1172
1173        int maxIdx = (int)data.mPolygons->size();
1174
1175        for (int i = 0; i < limit; ++ i)
1176        {
1177                // the already taken candidates are stored behind maxIdx
1178                // => assure that no index is taken twice
1179                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
1180                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
1181
1182                // swap candidate to the end to avoid testing same plane
1183                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
1184                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
1185
1186                // evaluate current candidate
1187                BspNodeGeometry fGeom, bGeom;
1188                float fArea, bArea;
1189                plane = poly->GetSupportingPlane();
1190                candidateCost = EvalSplitPlaneCost(plane, data, fGeom, bGeom, fArea, bArea);
1191               
1192                if (candidateCost < lowestCost)
1193                {
1194                        bestPlane = plane;
1195                        lowestCost = candidateCost;
1196                }
1197        }
1198
1199        // cost ratio miss
1200        if (mUsePolygonSplitIfAvailable && !data.mPolygons->empty())
1201        {
1202                frontData.mIsKdNode = backData.mIsKdNode = false;
1203                if (lowestCost > mTermMaxCostRatio)
1204                        return false;
1205
1206                return true;
1207        }
1208
1209        //-- evaluate axis aligned splits
1210        int axis;
1211        BspNodeGeometry *fGeom, *bGeom;
1212        float pFront, pBack;
1213
1214        candidateCost = SelectAxisAlignedPlane(plane,
1215                                                                                   data,
1216                                                                                   axis,
1217                                                                                   &fGeom,
1218                                                                                   &bGeom,
1219                                                                                   pFront,
1220                                                                                   pBack,
1221                                                                                   data.mIsKdNode);     
1222
1223        splitAxis = 3;
1224
1225        if (candidateCost < lowestCost)
1226        {       
1227                bestPlane = plane;
1228                lowestCost = candidateCost;
1229                splitAxis = axis;
1230                // assign already computed values
1231                // we can do this because we always save the
1232                // computed values from the axis aligned splits         
1233                if (fGeom && bGeom)
1234                {
1235                        frontData.mGeometry = fGeom;
1236                        backData.mGeometry = bGeom;
1237       
1238                        frontData.mProbability = pFront;
1239                        backData.mProbability = pBack;
1240                }
1241        }
1242        else
1243        {
1244                DEL_PTR(fGeom);
1245                DEL_PTR(bGeom);
1246        }
1247
1248        frontData.mIsKdNode = backData.mIsKdNode =
1249                (data.mIsKdNode && splitAxis < 3);
1250
1251#ifdef _DEBUG
1252        Debug << "plane lowest cost: " << lowestCost << endl;
1253#endif
1254
1255        // cost ratio miss
1256        if (lowestCost > mTermMaxCostRatio)
1257        {
1258                return false;
1259        }
1260
1261        return true;
1262}
1263
1264
1265Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
1266{
1267        const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1268
1269        const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
1270        const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();
1271
1272        const Vector3 pt = (maxPt + minPt) * 0.5;
1273        const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());
1274
1275        return Plane3(normal, pt);
1276}
1277
1278
1279Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
1280{
1281        Vector3 pt[3];
1282
1283        int idx[3];
1284        int cmaxT = 0;
1285        int cminT = 0;
1286        bool chooseMin = false;
1287
1288        for (int j = 0; j < 3; ++ j)
1289        {
1290                idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));
1291
1292                if (idx[j] >= (int)rays.size())
1293                {
1294                        idx[j] -= (int)rays.size();
1295
1296                        chooseMin = (cminT < 2);
1297                }
1298                else
1299                        chooseMin = (cmaxT < 2);
1300
1301                RayInfo rayInf = rays[idx[j]];
1302                pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
1303        }
1304
1305        return Plane3(pt[0], pt[1], pt[2]);
1306}
1307
1308
1309Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
1310{
1311        Vector3 pt[3];
1312
1313        int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1314        int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1315
1316        // check if rays different
1317        if (idx1 == idx2)
1318                idx2 = (idx2 + 1) % (int)rays.size();
1319
1320        const RayInfo ray1 = rays[idx1];
1321        const RayInfo ray2 = rays[idx2];
1322
1323        // normal vector of the plane parallel to both lines
1324        const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));
1325
1326        // vector from line 1 to line 2
1327        const Vector3 vd = ray2.ExtrapOrigin() - ray1.ExtrapOrigin();
1328
1329        // project vector on normal to get distance
1330        const float dist = DotProd(vd, norm);
1331
1332        // point on plane lies halfway between the two planes
1333        const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;
1334
1335        return Plane3(norm, planePt);
1336}
1337
1338
1339inline void VspBspTree::GenerateUniqueIdsForPvs()
1340{
1341        Intersectable::NewMail(); sBackId = Intersectable::sMailId;
1342        Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
1343        Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
1344}
1345
1346
1347float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
1348                                                                         const VspBspTraversalData &data,
1349                                                                         BspNodeGeometry &geomFront,
1350                                                                         BspNodeGeometry &geomBack,
1351                                                                         float &pFront,
1352                                                                         float &pBack) const
1353{
1354        float cost = 0;
1355
1356        float sumBalancedRays = 0;
1357        float sumRaySplits = 0;
1358
1359        int pvsFront = 0;
1360        int pvsBack = 0;
1361
1362        // probability that view point lies in back / front node
1363        float pOverall = 0;
1364        pFront = 0;
1365        pBack = 0;
1366
1367        int raysFront = 0;
1368        int raysBack = 0;
1369        int totalPvs = 0;
1370
1371        int limit;
1372        bool useRand;
1373
1374        // choose test rays randomly if too much
1375        if ((int)data.mRays->size() > mMaxTests)
1376        {
1377                useRand = true;
1378                limit = mMaxTests;
1379        }
1380        else
1381        {
1382                useRand = false;
1383                limit = (int)data.mRays->size();
1384        }
1385       
1386        for (int i = 0; i < limit; ++ i)
1387        {
1388                const int testIdx = useRand ?
1389                        (int)RandomValue(0, (Real)((int)data.mRays->size() - 1)) : i;
1390                RayInfo rayInf = (*data.mRays)[testIdx];
1391
1392                float t;
1393                VssRay *ray = rayInf.mRay;
1394                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
1395
1396        if (0)
1397                {
1398                        if (cf >= 0)
1399                                ++ raysFront;
1400                        if (cf <= 0)
1401                                ++ raysBack;
1402                }
1403
1404                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1405                {
1406                        sumBalancedRays += cf;
1407                }
1408
1409                if (mSplitPlaneStrategy & BALANCED_RAYS)
1410                {
1411                        if (cf == 0)
1412                                ++ sumRaySplits;
1413                }
1414
1415                if (mSplitPlaneStrategy & PVS)
1416                {
1417                        // find front and back pvs for origing and termination object
1418                        AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1419                        AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1420                }
1421        }
1422
1423        const float raysSize = (float)data.mRays->size() + Limits::Small;
1424
1425        if (mSplitPlaneStrategy & PVS)
1426        {
1427                // create unique ids for pvs heuristics
1428                GenerateUniqueIdsForPvs();
1429
1430                // construct child geometry with regard to the candidate split plane
1431                data.mGeometry->SplitGeometry(geomFront,
1432                                                                          geomBack,
1433                                                                          candidatePlane,
1434                                                                          mBox,
1435                                                                          mEpsilon);
1436
1437                pOverall = data.mProbability;
1438
1439                if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
1440                {
1441                        pFront = geomFront.GetVolume();
1442                        pBack = pOverall - pFront;
1443                }
1444                else
1445                {
1446                        pFront = geomFront.GetArea();
1447                        pBack = geomBack.GetArea();
1448                }
1449        }
1450
1451
1452        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1453                cost += mLeastRaySplitsFactor * sumRaySplits / raysSize;
1454
1455        if (mSplitPlaneStrategy & BALANCED_RAYS)
1456                cost += mBalancedRaysFactor * fabs(sumBalancedRays) / raysSize;
1457
1458        // -- pvs rendering heuristics
1459        if (mSplitPlaneStrategy & PVS)
1460        {
1461                const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1462                const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1463
1464                // only render cost heuristics or combined with standard deviation
1465                const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
1466        const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1467                const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1468                       
1469                const float oldRenderCost = pOverall * penaltyOld;
1470                const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1471
1472                float oldCost, newCost;
1473
1474                // only render cost
1475                if (1)
1476                {
1477                        oldCost = oldRenderCost;
1478                        newCost = newRenderCost;
1479                }
1480                else // also considering standard deviation
1481                {
1482                        // standard deviation is difference of back and front pvs
1483                        const float expectedCost = 0.5f * (penaltyFront + penaltyBack);
1484
1485                        const float newDeviation = 0.5f *
1486                                fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);
1487
1488                        const float oldDeviation = penaltyOld;
1489
1490                        newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
1491                        oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
1492                }
1493
1494                cost += mPvsFactor * newCost / (oldCost + Limits::Small);
1495               
1496        }
1497
1498#ifdef _DEBUG
1499        Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
1500                  << " frontpvs: " << pvsFront << " pFront: " << pFront
1501                  << " backpvs: " << pvsBack << " pBack: " << pBack << endl << endl;
1502        Debug << "cost: " << cost << endl;
1503#endif
1504
1505       
1506        // normalize cost by sum of linear factors
1507        if(0)
1508                return cost / (float)mCostNormalizer;
1509        else
1510                return cost;
1511}
1512
1513
1514float VspBspTree::EvalAxisAlignedSplitCost(const VspBspTraversalData &data,
1515                                                                                   const AxisAlignedBox3 &box,
1516                                                                                   const int axis,
1517                                                                                   const float &position,                                                                                 
1518                                                                                   float &pFront,
1519                                                                                   float &pBack) const
1520{
1521        int pvsTotal = 0;
1522        int pvsFront = 0;
1523        int pvsBack = 0;
1524       
1525        // create unique ids for pvs heuristics
1526        GenerateUniqueIdsForPvs();
1527
1528        const int pvsSize = data.mPvs;
1529
1530        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1531
1532        // this is the main ray classification loop!
1533        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
1534        {
1535                //if (!(*rit).mRay->IsActive()) continue;
1536
1537                // determine the side of this ray with respect to the plane
1538                float t;
1539                const int side = (*rit).ComputeRayIntersection(axis, position, t);
1540       
1541                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
1542                AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
1543        }
1544
1545        //-- pvs heuristics
1546        float pOverall;
1547
1548        //-- compute heurstics
1549        //   we take simplified computation for mid split
1550               
1551        pOverall = data.mProbability;
1552
1553        if (!mUseAreaForPvs)
1554        {   // volume
1555                pBack = pFront = pOverall * 0.5f;
1556               
1557#if 0
1558                // box length substitute for probability
1559                const float minBox = box.Min(axis);
1560                const float maxBox = box.Max(axis);
1561
1562                pBack = position - minBox;
1563                pFront = maxBox - position;
1564                pOverall = maxBox - minBox;
1565#endif
1566        }
1567        else //-- area substitute for probability
1568        {
1569                const int axis2 = (axis + 1) % 3;
1570                const int axis3 = (axis + 2) % 3;
1571
1572                const float faceArea =
1573                        (box.Max(axis2) - box.Min(axis2)) *
1574                        (box.Max(axis3) - box.Min(axis3));
1575
1576                pBack = pFront = pOverall * 0.5f + faceArea;
1577        }
1578
1579#ifdef _DEBUG
1580        Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
1581        Debug << pFront << " " << pBack << " " << pOverall << endl;
1582#endif
1583
1584       
1585        const float newCost = pvsBack * pBack + pvsFront * pFront;
1586        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
1587
1588        return  (mCtDivCi + newCost) / oldCost;
1589}
1590
1591
1592void VspBspTree::AddObjToPvs(Intersectable *obj,
1593                                                         const int cf,
1594                                                         int &frontPvs,
1595                                                         int &backPvs,
1596                                                         int &totalPvs) const
1597{
1598        if (!obj)
1599                return;
1600       
1601        if ((obj->mMailbox != sFrontId) &&
1602                (obj->mMailbox != sBackId) &&
1603                (obj->mMailbox != sFrontAndBackId))
1604        {
1605                ++ totalPvs;
1606        }
1607
1608        // TODO: does this really belong to no pvs?
1609        //if (cf == Ray::COINCIDENT) return;
1610
1611        // object belongs to both PVS
1612        if (cf >= 0)
1613        {
1614                if ((obj->mMailbox != sFrontId) &&
1615                        (obj->mMailbox != sFrontAndBackId))
1616                {
1617                        ++ frontPvs;
1618               
1619                        if (obj->mMailbox == sBackId)
1620                                obj->mMailbox = sFrontAndBackId;
1621                        else
1622                                obj->mMailbox = sFrontId;
1623                }
1624        }
1625
1626        if (cf <= 0)
1627        {
1628                if ((obj->mMailbox != sBackId) &&
1629                        (obj->mMailbox != sFrontAndBackId))
1630                {
1631                        ++ backPvs;
1632               
1633                        if (obj->mMailbox == sFrontId)
1634                                obj->mMailbox = sFrontAndBackId;
1635                        else
1636                                obj->mMailbox = sBackId;
1637                }
1638        }
1639}
1640
1641
1642void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves,
1643                                                           const bool onlyUnmailed,
1644                                                           const int maxPvsSize) const
1645{
1646        stack<BspNode *> nodeStack;
1647        nodeStack.push(mRoot);
1648
1649        while (!nodeStack.empty())
1650        {
1651                BspNode *node = nodeStack.top();
1652                nodeStack.pop();
1653               
1654                if (node->IsLeaf())
1655                {
1656                        // test if this leaf is in valid view space
1657                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1658                        if (leaf->TreeValid() &&
1659                                (!onlyUnmailed || !leaf->Mailed()) &&
1660                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
1661                        {
1662                                leaves.push_back(leaf);
1663                        }
1664                }
1665                else
1666                {
1667                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1668
1669                        nodeStack.push(interior->GetBack());
1670                        nodeStack.push(interior->GetFront());
1671                }
1672        }
1673}
1674
1675
1676AxisAlignedBox3 VspBspTree::GetBoundingBox() const
1677{
1678        return mBox;
1679}
1680
1681
1682BspNode *VspBspTree::GetRoot() const
1683{
1684        return mRoot;
1685}
1686
1687
1688void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
1689{
1690        // the node became a leaf -> evaluate stats for leafs
1691        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
1692
1693        // store maximal and minimal depth
1694        if (data.mDepth > mBspStats.maxDepth)
1695                mBspStats.maxDepth = data.mDepth;
1696
1697        if (data.mPvs > mBspStats.maxPvs)
1698                mBspStats.maxPvs = data.mPvs;
1699       
1700        if (data.mDepth < mBspStats.minDepth)
1701                mBspStats.minDepth = data.mDepth;
1702
1703        if (data.mDepth >= mTermMaxDepth)
1704                ++ mBspStats.maxDepthNodes;
1705
1706        // accumulate rays to compute rays /  leaf
1707        mBspStats.accumRays += (int)data.mRays->size();
1708
1709        if (data.mPvs < mTermMinPvs)
1710                ++ mBspStats.minPvsNodes;
1711
1712        if ((int)data.mRays->size() < mTermMinRays)
1713                ++ mBspStats.minRaysNodes;
1714
1715        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1716                ++ mBspStats.maxRayContribNodes;
1717
1718        if (data.mProbability <= mTermMinProbability)
1719                ++ mBspStats.minProbabilityNodes;
1720       
1721        // accumulate depth to compute average depth
1722        mBspStats.accumDepth += data.mDepth;
1723
1724        ++ mCreatedViewCells;
1725#ifdef _DEBUG
1726        Debug << "BSP stats: "
1727                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1728                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1729          //              << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
1730                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1731                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
1732                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1733#endif
1734}
1735
1736
1737int VspBspTree::CastRay(Ray &ray)
1738{
1739        int hits = 0;
1740
1741        stack<BspRayTraversalData> tQueue;
1742
1743        float maxt, mint;
1744
1745        if (!mBox.GetRaySegment(ray, mint, maxt))
1746                return 0;
1747
1748        Intersectable::NewMail();
1749        ViewCell::NewMail();
1750        Vector3 entp = ray.Extrap(mint);
1751        Vector3 extp = ray.Extrap(maxt);
1752
1753        BspNode *node = mRoot;
1754        BspNode *farChild = NULL;
1755
1756        while (1)
1757        {
1758                if (!node->IsLeaf())
1759                {
1760                        BspInterior *in = dynamic_cast<BspInterior *>(node);
1761
1762                        Plane3 splitPlane = in->GetPlane();
1763                        const int entSide = splitPlane.Side(entp);
1764                        const int extSide = splitPlane.Side(extp);
1765
1766                        if (entSide < 0)
1767                        {
1768                                node = in->GetBack();
1769
1770                                if(extSide <= 0) // plane does not split ray => no far child
1771                                        continue;
1772
1773                                farChild = in->GetFront(); // plane splits ray
1774
1775                        } else if (entSide > 0)
1776                        {
1777                                node = in->GetFront();
1778
1779                                if (extSide >= 0) // plane does not split ray => no far child
1780                                        continue;
1781
1782                                farChild = in->GetBack(); // plane splits ray
1783                        }
1784                        else // ray and plane are coincident
1785                        {
1786                                // WHAT TO DO IN THIS CASE ?
1787                                //break;
1788                                node = in->GetFront();
1789                                continue;
1790                        }
1791
1792                        // push data for far child
1793                        tQueue.push(BspRayTraversalData(farChild, extp, maxt));
1794
1795                        // find intersection of ray segment with plane
1796                        float t;
1797                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
1798                        maxt *= t;
1799
1800                } else // reached leaf => intersection with view cell
1801                {
1802                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1803
1804                        if (!leaf->GetViewCell()->Mailed())
1805                        {
1806                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
1807                                leaf->GetViewCell()->Mail();
1808                                ++ hits;
1809                        }
1810
1811                        //-- fetch the next far child from the stack
1812                        if (tQueue.empty())
1813                                break;
1814
1815                        entp = extp;
1816                        mint = maxt; // NOTE: need this?
1817
1818                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
1819                                break;
1820
1821                        BspRayTraversalData &s = tQueue.top();
1822
1823                        node = s.mNode;
1824                        extp = s.mExitPoint;
1825                        maxt = s.mMaxT;
1826
1827                        tQueue.pop();
1828                }
1829        }
1830
1831        return hits;
1832}
1833
1834
1835void VspBspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
1836{
1837        ViewCell::NewMail();
1838       
1839        CollectViewCells(mRoot, onlyValid, viewCells, true);
1840}
1841
1842
1843void VspBspTree::CollapseViewCells()
1844{
1845// TODO
1846#if VC_HISTORY
1847        stack<BspNode *> nodeStack;
1848
1849        if (!mRoot)
1850                return;
1851
1852        nodeStack.push(mRoot);
1853       
1854        while (!nodeStack.empty())
1855        {
1856                BspNode *node = nodeStack.top();
1857                nodeStack.pop();
1858               
1859                if (node->IsLeaf())
1860        {
1861                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
1862
1863                        if (!viewCell->GetValid())
1864                        {
1865                                BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
1866       
1867                                ViewCellContainer leaves;
1868                                mViewCellsTree->CollectLeaves(viewCell, leaves);
1869
1870                                ViewCellContainer::const_iterator it, it_end = leaves.end();
1871
1872                                for (it = leaves.begin(); it != it_end; ++ it)
1873                                {
1874                                        BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
1875                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
1876                                        ++ mBspStats.invalidLeaves;
1877                                }
1878
1879                                // add to unbounded view cell
1880                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
1881                                DEL_PTR(viewCell);
1882                        }
1883                }
1884                else
1885                {
1886                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1887               
1888                        nodeStack.push(interior->GetFront());
1889                        nodeStack.push(interior->GetBack());
1890                }
1891        }
1892
1893        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
1894#endif
1895}
1896
1897
1898void VspBspTree::CollectRays(VssRayContainer &rays)
1899{
1900        vector<BspLeaf *> leaves;
1901
1902        vector<BspLeaf *>::const_iterator lit, lit_end = leaves.end();
1903
1904        for (lit = leaves.begin(); lit != lit_end; ++ lit)
1905        {
1906                BspLeaf *leaf = *lit;
1907                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
1908
1909                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
1910                        rays.push_back(*rit);
1911        }
1912}
1913
1914
1915void VspBspTree::ValidateTree()
1916{
1917        stack<BspNode *> nodeStack;
1918
1919        if (!mRoot)
1920                return;
1921
1922        nodeStack.push(mRoot);
1923       
1924        mBspStats.invalidLeaves = 0;
1925        while (!nodeStack.empty())
1926        {
1927                BspNode *node = nodeStack.top();
1928                nodeStack.pop();
1929               
1930                if (node->IsLeaf())
1931                {
1932                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1933
1934                        if (!leaf->GetViewCell()->GetValid())
1935                                ++ mBspStats.invalidLeaves;
1936
1937                        // validity flags don't match => repair
1938                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
1939                        {
1940                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
1941                                PropagateUpValidity(leaf);
1942                        }
1943                }
1944                else
1945                {
1946                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1947               
1948                        nodeStack.push(interior->GetFront());
1949                        nodeStack.push(interior->GetBack());
1950                }
1951        }
1952
1953        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
1954}
1955
1956
1957
1958void VspBspTree::CollectViewCells(BspNode *root,
1959                                                                  bool onlyValid,
1960                                                                  ViewCellContainer &viewCells,
1961                                                                  bool onlyUnmailed) const
1962{
1963        stack<BspNode *> nodeStack;
1964
1965        if (!root)
1966                return;
1967
1968        nodeStack.push(root);
1969       
1970        while (!nodeStack.empty())
1971        {
1972                BspNode *node = nodeStack.top();
1973                nodeStack.pop();
1974               
1975                if (node->IsLeaf())
1976                {
1977                        if (!onlyValid || node->TreeValid())
1978                        {
1979                                ViewCell *viewCell =
1980                                        mViewCellsTree->GetActiveViewCell(dynamic_cast<BspLeaf *>(node)->GetViewCell());
1981                                               
1982                                if (!onlyUnmailed || !viewCell->Mailed())
1983                                {
1984                                        viewCell->Mail();
1985                                        viewCells.push_back(viewCell);
1986                                }
1987                        }
1988                }
1989                else
1990                {
1991                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1992               
1993                        nodeStack.push(interior->GetFront());
1994                        nodeStack.push(interior->GetBack());
1995                }
1996        }
1997
1998}
1999
2000
2001void VspBspTree::PreprocessPolygons(PolygonContainer &polys)
2002{
2003        // preprocess: throw out polygons coincident to the view space box (not needed)
2004        PolygonContainer boxPolys;
2005        mBox.ExtractPolys(boxPolys);
2006        vector<Plane3> boxPlanes;
2007
2008        PolygonContainer::iterator pit, pit_end = boxPolys.end();
2009
2010        // extract planes of box
2011        // TODO: can be done more elegantly than first extracting polygons
2012        // and take their planes
2013        for (pit = boxPolys.begin(); pit != pit_end; ++ pit)
2014        {
2015                boxPlanes.push_back((*pit)->GetSupportingPlane());
2016        }
2017
2018        pit_end = polys.end();
2019
2020        for (pit = polys.begin(); pit != pit_end; ++ pit)
2021        {
2022                vector<Plane3>::const_iterator bit, bit_end = boxPlanes.end();
2023               
2024                for (bit = boxPlanes.begin(); (bit != bit_end) && (*pit); ++ bit)
2025                {
2026                        const int cf = (*pit)->ClassifyPlane(*bit, mEpsilon);
2027
2028                        if (cf == Polygon3::COINCIDENT)
2029                        {
2030                                DEL_PTR(*pit);
2031                                //Debug << "coincident!!" << endl;
2032                        }
2033                }
2034        }
2035
2036        // remove deleted entries
2037        for (int i = 0; i < (int)polys.size(); ++ i)
2038        {
2039                while (!polys[i] && (i < (int)polys.size()))
2040                {
2041                        swap(polys[i], polys.back());
2042                        polys.pop_back();
2043                }
2044        }
2045}
2046
2047
2048float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
2049{
2050        float len = 0;
2051
2052        RayInfoContainer::const_iterator it, it_end = rays.end();
2053
2054        for (it = rays.begin(); it != it_end; ++ it)
2055                len += (*it).SegmentLength();
2056
2057        return len;
2058}
2059
2060
2061int VspBspTree::SplitRays(const Plane3 &plane,
2062                                                  RayInfoContainer &rays,
2063                                                  RayInfoContainer &frontRays,
2064                                                  RayInfoContainer &backRays) const
2065{
2066        int splits = 0;
2067
2068        RayInfoContainer::const_iterator it, it_end = rays.end();
2069
2070        for (it = rays.begin(); it != it_end; ++ it)
2071        {
2072                RayInfo bRay = *it;
2073               
2074                VssRay *ray = bRay.mRay;
2075                float t;
2076
2077                // get classification and receive new t
2078                const int cf = bRay.ComputeRayIntersection(plane, t);
2079
2080                switch (cf)
2081                {
2082                case -1:
2083                        backRays.push_back(bRay);
2084                        break;
2085                case 1:
2086                        frontRays.push_back(bRay);
2087                        break;
2088                case 0:
2089                        {
2090                                //-- split ray
2091                                //   test if start point behind or in front of plane
2092                                const int side = plane.Side(bRay.ExtrapOrigin());
2093
2094                                if (side <= 0)
2095                                {
2096                                        backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2097                                        frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2098                                }
2099                                else
2100                                {
2101                                        frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2102                                        backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2103                                }
2104                        }
2105                        break;
2106                default:
2107                        Debug << "Should not come here" << endl;
2108                        break;
2109                }
2110        }
2111
2112        return splits;
2113}
2114
2115
2116void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
2117{
2118        BspNode *lastNode;
2119
2120        do
2121        {
2122                lastNode = n;
2123
2124                // want to get planes defining geometry of this node => don't take
2125                // split plane of node itself
2126                n = n->GetParent();
2127
2128                if (n)
2129                {
2130                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
2131                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
2132
2133            if (interior->GetFront() != lastNode)
2134                                halfSpace.ReverseOrientation();
2135
2136                        halfSpaces.push_back(halfSpace);
2137                }
2138        }
2139        while (n);
2140}
2141
2142
2143void VspBspTree::ConstructGeometry(BspNode *n,
2144                                                                   BspNodeGeometry &geom) const
2145{
2146        vector<Plane3> halfSpaces;
2147        ExtractHalfSpaces(n, halfSpaces);
2148
2149        PolygonContainer candidatePolys;
2150
2151        // bounded planes are added to the polygons (reverse polygons
2152        // as they have to be outfacing
2153        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
2154        {
2155                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
2156
2157                if (p->Valid(mEpsilon))
2158                {
2159                        candidatePolys.push_back(p->CreateReversePolygon());
2160                        DEL_PTR(p);
2161                }
2162        }
2163
2164        // add faces of bounding box (also could be faces of the cell)
2165        for (int i = 0; i < 6; ++ i)
2166        {
2167                VertexContainer vertices;
2168
2169                for (int j = 0; j < 4; ++ j)
2170                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
2171
2172                candidatePolys.push_back(new Polygon3(vertices));
2173        }
2174
2175        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
2176        {
2177                // polygon is split by all other planes
2178                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
2179                {
2180                        if (i == j) // polygon and plane are coincident
2181                                continue;
2182
2183                        VertexContainer splitPts;
2184                        Polygon3 *frontPoly, *backPoly;
2185
2186                        const int cf =
2187                                candidatePolys[i]->ClassifyPlane(halfSpaces[j],
2188                                                                                                 mEpsilon);
2189
2190                        switch (cf)
2191                        {
2192                                case Polygon3::SPLIT:
2193                                        frontPoly = new Polygon3();
2194                                        backPoly = new Polygon3();
2195
2196                                        candidatePolys[i]->Split(halfSpaces[j],
2197                                                                                         *frontPoly,
2198                                                                                         *backPoly,
2199                                                                                         mEpsilon);
2200
2201                                        DEL_PTR(candidatePolys[i]);
2202
2203                                        if (frontPoly->Valid(mEpsilon))
2204                                                candidatePolys[i] = frontPoly;
2205                                        else
2206                                                DEL_PTR(frontPoly);
2207
2208                                        DEL_PTR(backPoly);
2209                                        break;
2210                                case Polygon3::BACK_SIDE:
2211                                        DEL_PTR(candidatePolys[i]);
2212                                        break;
2213                                // just take polygon as it is
2214                                case Polygon3::FRONT_SIDE:
2215                                case Polygon3::COINCIDENT:
2216                                default:
2217                                        break;
2218                        }
2219                }
2220
2221                if (candidatePolys[i])
2222                        geom.mPolys.push_back(candidatePolys[i]);
2223        }
2224}
2225
2226
2227void VspBspTree::ConstructGeometry(ViewCell *vc,
2228                                                                   BspNodeGeometry &vcGeom) const
2229{
2230        ViewCellContainer leaves;
2231       
2232        mViewCellsTree->CollectLeaves(vc, leaves);
2233
2234        ViewCellContainer::const_iterator it, it_end = leaves.end();
2235
2236        for (it = leaves.begin(); it != it_end; ++ it)
2237        {
2238                BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
2239               
2240                ConstructGeometry(l, vcGeom);
2241        }
2242}
2243
2244
2245typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;
2246
2247
2248int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
2249                                                          const bool onlyUnmailed) const
2250{
2251        stack<bspNodePair> nodeStack;
2252       
2253        BspNodeGeometry nodeGeom;
2254        ConstructGeometry(n, nodeGeom);
2255       
2256        // split planes from the root to this node
2257        // needed to verify that we found neighbor leaf
2258        // TODO: really needed?
2259        vector<Plane3> halfSpaces;
2260        ExtractHalfSpaces(n, halfSpaces);
2261
2262
2263        BspNodeGeometry *rgeom = new BspNodeGeometry();
2264        ConstructGeometry(mRoot, *rgeom);
2265
2266        nodeStack.push(bspNodePair(mRoot, rgeom));
2267
2268        while (!nodeStack.empty())
2269        {
2270                BspNode *node = nodeStack.top().first;
2271                BspNodeGeometry *geom = nodeStack.top().second;
2272       
2273                nodeStack.pop();
2274
2275                if (node->IsLeaf())
2276                {
2277                        // test if this leaf is in valid view space
2278                        if (node->TreeValid() &&
2279                                (node != n) &&
2280                                (!onlyUnmailed || !node->Mailed()))
2281                        {
2282                                bool isAdjacent = true;
2283
2284                                if (1)
2285                                {
2286                                        // test all planes of current node if still adjacent
2287                                        for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
2288                                        {
2289                                                const int cf =
2290                                                        Polygon3::ClassifyPlane(geom->mPolys,
2291                                                                                                        halfSpaces[i],
2292                                                                                                        mEpsilon);
2293
2294                                                if (cf == Polygon3::BACK_SIDE)
2295                                                {
2296                                                        isAdjacent = false;
2297                                                }
2298                                        }
2299                                }
2300                                else
2301                                {
2302                                        // TODO: why is this wrong??
2303                                        // test all planes of current node if still adjacent
2304                                        for (int i = 0; (i < (int)nodeGeom.mPolys.size()) && isAdjacent; ++ i)
2305                                        {
2306                                                Polygon3 *poly = nodeGeom.mPolys[i];
2307
2308                                                const int cf =
2309                                                        Polygon3::ClassifyPlane(geom->mPolys,
2310                                                                                                        poly->GetSupportingPlane(),
2311                                                                                                        mEpsilon);
2312
2313                                                if (cf == Polygon3::BACK_SIDE)
2314                                                {
2315                                                        isAdjacent = false;
2316                                                }
2317                                        }
2318                                }
2319                                // neighbor was found
2320                                if (isAdjacent)
2321                                {       
2322                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2323                                }
2324                        }
2325                }
2326                else
2327                {
2328                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2329
2330                        const int cf = Polygon3::ClassifyPlane(nodeGeom.mPolys,
2331                                                                                                   interior->GetPlane(),
2332                                                                                                   mEpsilon);
2333                       
2334                        BspNode *front = interior->GetFront();
2335                        BspNode *back = interior->GetBack();
2336           
2337                        BspNodeGeometry *fGeom = new BspNodeGeometry();
2338                        BspNodeGeometry *bGeom = new BspNodeGeometry();
2339
2340                        geom->SplitGeometry(*fGeom,
2341                                                                *bGeom,
2342                                                                interior->GetPlane(),
2343                                                                mBox,
2344                                                                mEpsilon);
2345               
2346                        if (cf == Polygon3::FRONT_SIDE)
2347                        {
2348                                nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2349                                DEL_PTR(bGeom);
2350                        }
2351                        else
2352                        {
2353                                if (cf == Polygon3::BACK_SIDE)
2354                                {
2355                                        nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2356                                        DEL_PTR(fGeom);
2357                                }
2358                                else
2359                                {       // random decision
2360                                        nodeStack.push(bspNodePair(front, fGeom));
2361                                        nodeStack.push(bspNodePair(back, bGeom));
2362                                }
2363                        }
2364                }
2365       
2366                DEL_PTR(geom);
2367        }
2368
2369        return (int)neighbors.size();
2370}
2371
2372
2373
2374int VspBspTree::FindApproximateNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
2375                                                                                 const bool onlyUnmailed) const
2376{
2377        stack<bspNodePair> nodeStack;
2378       
2379        BspNodeGeometry nodeGeom;
2380        ConstructGeometry(n, nodeGeom);
2381       
2382        // split planes from the root to this node
2383        // needed to verify that we found neighbor leaf
2384        // TODO: really needed?
2385        vector<Plane3> halfSpaces;
2386        ExtractHalfSpaces(n, halfSpaces);
2387
2388
2389        BspNodeGeometry *rgeom = new BspNodeGeometry();
2390        ConstructGeometry(mRoot, *rgeom);
2391
2392        nodeStack.push(bspNodePair(mRoot, rgeom));
2393
2394        while (!nodeStack.empty())
2395        {
2396                BspNode *node = nodeStack.top().first;
2397                BspNodeGeometry *geom = nodeStack.top().second;
2398       
2399                nodeStack.pop();
2400
2401                if (node->IsLeaf())
2402                {
2403                        // test if this leaf is in valid view space
2404                        if (node->TreeValid() &&
2405                                (node != n) &&
2406                                (!onlyUnmailed || !node->Mailed()))
2407                        {
2408                                bool isAdjacent = true;
2409
2410                                // neighbor was found
2411                                if (isAdjacent)
2412                                {       
2413                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2414                                }
2415                        }
2416                }
2417                else
2418                {
2419                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2420
2421                        const int cf = Polygon3::ClassifyPlane(nodeGeom.mPolys,
2422                                                                                                   interior->GetPlane(),
2423                                                                                                   mEpsilon);
2424                       
2425                        BspNode *front = interior->GetFront();
2426                        BspNode *back = interior->GetBack();
2427           
2428                        BspNodeGeometry *fGeom = new BspNodeGeometry();
2429                        BspNodeGeometry *bGeom = new BspNodeGeometry();
2430
2431                        geom->SplitGeometry(*fGeom,
2432                                                                *bGeom,
2433                                                                interior->GetPlane(),
2434                                                                mBox,
2435                                                                mEpsilon);
2436               
2437                        if (cf == Polygon3::FRONT_SIDE)
2438                        {
2439                                nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2440                                DEL_PTR(bGeom);
2441                        }
2442                        else
2443                        {
2444                                if (cf == Polygon3::BACK_SIDE)
2445                                {
2446                                        nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2447                                        DEL_PTR(fGeom);
2448                                }
2449                                else
2450                                {       // random decision
2451                                        nodeStack.push(bspNodePair(front, fGeom));
2452                                        nodeStack.push(bspNodePair(back, bGeom));
2453                                }
2454                        }
2455                }
2456       
2457                DEL_PTR(geom);
2458        }
2459
2460        return (int)neighbors.size();
2461}
2462
2463
2464
2465BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
2466{
2467    stack<BspNode *> nodeStack;
2468        nodeStack.push(mRoot);
2469
2470        int mask = rand();
2471
2472        while (!nodeStack.empty())
2473        {
2474                BspNode *node = nodeStack.top();
2475                nodeStack.pop();
2476
2477                if (node->IsLeaf())
2478                {
2479                        return dynamic_cast<BspLeaf *>(node);
2480                }
2481                else
2482                {
2483                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2484                        BspNode *next;
2485                        BspNodeGeometry geom;
2486
2487                        // todo: not very efficient: constructs full cell everytime
2488                        ConstructGeometry(interior, geom);
2489
2490                        const int cf =
2491                                Polygon3::ClassifyPlane(geom.mPolys, halfspace, mEpsilon);
2492
2493                        if (cf == Polygon3::BACK_SIDE)
2494                                next = interior->GetFront();
2495                        else
2496                                if (cf == Polygon3::FRONT_SIDE)
2497                                        next = interior->GetFront();
2498                        else
2499                        {
2500                                // random decision
2501                                if (mask & 1)
2502                                        next = interior->GetBack();
2503                                else
2504                                        next = interior->GetFront();
2505                                mask = mask >> 1;
2506                        }
2507
2508                        nodeStack.push(next);
2509                }
2510        }
2511
2512        return NULL;
2513}
2514
2515BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
2516{
2517        stack<BspNode *> nodeStack;
2518
2519        nodeStack.push(mRoot);
2520
2521        int mask = rand();
2522
2523        while (!nodeStack.empty())
2524        {
2525                BspNode *node = nodeStack.top();
2526                nodeStack.pop();
2527
2528                if (node->IsLeaf())
2529                {
2530                        if ( (!onlyUnmailed || !node->Mailed()) )
2531                                return dynamic_cast<BspLeaf *>(node);
2532                }
2533                else
2534                {
2535                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2536
2537                        // random decision
2538                        if (mask & 1)
2539                                nodeStack.push(interior->GetBack());
2540                        else
2541                                nodeStack.push(interior->GetFront());
2542
2543                        mask = mask >> 1;
2544                }
2545        }
2546
2547        return NULL;
2548}
2549
2550int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
2551{
2552        int pvsSize = 0;
2553
2554        RayInfoContainer::const_iterator rit, rit_end = rays.end();
2555
2556        Intersectable::NewMail();
2557
2558        for (rit = rays.begin(); rit != rays.end(); ++ rit)
2559        {
2560                VssRay *ray = (*rit).mRay;
2561
2562                if (ray->mOriginObject)
2563                {
2564                        if (!ray->mOriginObject->Mailed())
2565                        {
2566                                ray->mOriginObject->Mail();
2567                                ++ pvsSize;
2568                        }
2569                }
2570                if (ray->mTerminationObject)
2571                {
2572                        if (!ray->mTerminationObject->Mailed())
2573                        {
2574                                ray->mTerminationObject->Mail();
2575                                ++ pvsSize;
2576                        }
2577                }
2578        }
2579
2580        return pvsSize;
2581}
2582
2583float VspBspTree::GetEpsilon() const
2584{
2585        return mEpsilon;
2586}
2587
2588
2589int VspBspTree::SplitPolygons(const Plane3 &plane,
2590                                                          PolygonContainer &polys,
2591                                                          PolygonContainer &frontPolys,
2592                                                          PolygonContainer &backPolys,
2593                                                          PolygonContainer &coincident) const
2594{
2595        int splits = 0;
2596
2597        PolygonContainer::const_iterator it, it_end = polys.end();
2598
2599        for (it = polys.begin(); it != polys.end(); ++ it)     
2600        {
2601                Polygon3 *poly = *it;
2602
2603                // classify polygon
2604                const int cf = poly->ClassifyPlane(plane, mEpsilon);
2605
2606                switch (cf)
2607                {
2608                        case Polygon3::COINCIDENT:
2609                                coincident.push_back(poly);
2610                                break;
2611                        case Polygon3::FRONT_SIDE:
2612                                frontPolys.push_back(poly);
2613                                break;
2614                        case Polygon3::BACK_SIDE:
2615                                backPolys.push_back(poly);
2616                                break;
2617                        case Polygon3::SPLIT:
2618                                backPolys.push_back(poly);
2619                                frontPolys.push_back(poly);
2620                                ++ splits;
2621                                break;
2622                        default:
2623                Debug << "SHOULD NEVER COME HERE\n";
2624                                break;
2625                }
2626        }
2627
2628        return splits;
2629}
2630
2631
2632int VspBspTree::CastLineSegment(const Vector3 &origin,
2633                                                                const Vector3 &termination,
2634                                                                vector<ViewCell *> &viewcells)
2635{
2636        int hits = 0;
2637        stack<BspRayTraversalData> tQueue;
2638
2639        float mint = 0.0f, maxt = 1.0f;
2640
2641        Intersectable::NewMail();
2642        ViewCell::NewMail();
2643
2644        Vector3 entp = origin;
2645        Vector3 extp = termination;
2646
2647        BspNode *node = mRoot;
2648        BspNode *farChild = NULL;
2649
2650        float t;
2651        while (1)
2652        {
2653                if (!node->IsLeaf())
2654                {
2655                        BspInterior *in = dynamic_cast<BspInterior *>(node);
2656
2657                        Plane3 splitPlane = in->GetPlane();
2658                       
2659                        const int entSide = splitPlane.Side(entp);
2660                        const int extSide = splitPlane.Side(extp);
2661
2662                        if (entSide < 0)
2663                        {
2664                          node = in->GetBack();
2665                          // plane does not split ray => no far child
2666                          if (extSide <= 0)
2667                                continue;
2668                         
2669                          farChild = in->GetFront(); // plane splits ray
2670                        }
2671                        else if (entSide > 0)
2672                        {
2673                          node = in->GetFront();
2674
2675                          if (extSide >= 0) // plane does not split ray => no far child
2676                                continue;
2677
2678                                farChild = in->GetBack(); // plane splits ray
2679                        }
2680                        else // ray end point on plane
2681                        {       // NOTE: what to do if ray is coincident with plane?
2682                                if (extSide < 0)
2683                                        node = in->GetBack();
2684                                else
2685                                        node = in->GetFront();
2686                                                               
2687                                continue; // no far child
2688                        }
2689
2690                        // push data for far child
2691                        tQueue.push(BspRayTraversalData(farChild, extp));
2692
2693                        // find intersection of ray segment with plane
2694                        extp = splitPlane.FindIntersection(origin, extp, &t);
2695                }
2696                else
2697                {
2698                        // reached leaf => intersection with view cell
2699                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2700
2701                        ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
2702                        if (!viewCell->Mailed())
2703                        {
2704                                viewcells.push_back(viewCell);
2705                                viewCell->Mail();
2706                                ++ hits;
2707                        }
2708
2709                        //-- fetch the next far child from the stack
2710                        if (tQueue.empty())
2711                                break;
2712
2713                        entp = extp;
2714                       
2715                        BspRayTraversalData &s = tQueue.top();
2716
2717                        node = s.mNode;
2718                        extp = s.mExitPoint;
2719
2720                        tQueue.pop();
2721                }
2722        }
2723
2724        return hits;
2725}
2726
2727
2728
2729
2730int VspBspTree::TreeDistance(BspNode *n1, BspNode *n2) const
2731{
2732        std::deque<BspNode *> path1;
2733        BspNode *p1 = n1;
2734
2735        // create path from node 1 to root
2736        while (p1)
2737        {
2738                if (p1 == n2) // second node on path
2739                        return (int)path1.size();
2740
2741                path1.push_front(p1);
2742                p1 = p1->GetParent();
2743        }
2744
2745        int depth = n2->GetDepth();
2746        int d = depth;
2747
2748        BspNode *p2 = n2;
2749
2750        // compare with same depth
2751        while (1)
2752        {
2753                if ((d < (int)path1.size()) && (p2 == path1[d]))
2754                        return (depth - d) + ((int)path1.size() - 1 - d);
2755
2756                -- d;
2757                p2 = p2->GetParent();
2758        }
2759
2760        return 0; // never come here
2761}
2762
2763
2764BspNode *VspBspTree::CollapseTree(BspNode *node, int &collapsed)
2765{
2766// TODO
2767#if VC_HISTORY
2768        if (node->IsLeaf())
2769                return node;
2770
2771        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2772
2773        BspNode *front = CollapseTree(interior->GetFront(), collapsed);
2774        BspNode *back = CollapseTree(interior->GetBack(), collapsed);
2775
2776        if (front->IsLeaf() && back->IsLeaf())
2777        {
2778                BspLeaf *frontLeaf = dynamic_cast<BspLeaf *>(front);
2779                BspLeaf *backLeaf = dynamic_cast<BspLeaf *>(back);
2780
2781                //-- collapse tree
2782                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
2783                {
2784                        BspViewCell *vc = frontLeaf->GetViewCell();
2785
2786                        BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);
2787                        leaf->SetTreeValid(frontLeaf->TreeValid());
2788
2789                        // replace a link from node's parent
2790                        if (leaf->GetParent())
2791                                leaf->GetParent()->ReplaceChildLink(node, leaf);
2792                        else
2793                                mRoot = leaf;
2794
2795                        ++ collapsed;
2796                        delete interior;
2797
2798                        return leaf;
2799                }
2800        }
2801#endif
2802        return node;
2803}
2804
2805
2806int VspBspTree::CollapseTree()
2807{
2808        int collapsed = 0;
2809        //TODO
2810#if VC_HISTORY
2811        (void)CollapseTree(mRoot, collapsed);
2812
2813        // revalidate leaves
2814        RepairViewCellsLeafLists();
2815#endif
2816        return collapsed;
2817}
2818
2819
2820void VspBspTree::RepairViewCellsLeafLists()
2821{
2822// TODO
2823#if VC_HISTORY
2824        // list not valid anymore => clear
2825        stack<BspNode *> nodeStack;
2826        nodeStack.push(mRoot);
2827
2828        ViewCell::NewMail();
2829
2830        while (!nodeStack.empty())
2831        {
2832                BspNode *node = nodeStack.top();
2833                nodeStack.pop();
2834
2835                if (node->IsLeaf())
2836                {
2837                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2838
2839                        BspViewCell *viewCell = leaf->GetViewCell();
2840
2841                        if (!viewCell->Mailed())
2842                        {
2843                                viewCell->mLeaves.clear();
2844                                viewCell->Mail();
2845                        }
2846       
2847                        viewCell->mLeaves.push_back(leaf);
2848
2849                }
2850                else
2851                {
2852                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2853
2854                        nodeStack.push(interior->GetFront());
2855                        nodeStack.push(interior->GetBack());
2856                }
2857        }
2858// TODO
2859#endif
2860}
2861
2862
2863typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;
2864
2865
2866int VspBspTree::CastBeam(Beam &beam)
2867{
2868    stack<bspNodePair> nodeStack;
2869        BspNodeGeometry *rgeom = new BspNodeGeometry();
2870        ConstructGeometry(mRoot, *rgeom);
2871
2872        nodeStack.push(bspNodePair(mRoot, rgeom));
2873 
2874        ViewCell::NewMail();
2875
2876        while (!nodeStack.empty())
2877        {
2878                BspNode *node = nodeStack.top().first;
2879                BspNodeGeometry *geom = nodeStack.top().second;
2880                nodeStack.pop();
2881               
2882                AxisAlignedBox3 box;
2883                box.Initialize();
2884                geom->IncludeInBox(box);
2885
2886                const int side = beam.ComputeIntersection(box);
2887               
2888                switch (side)
2889                {
2890                case -1:
2891                        CollectViewCells(node, true, beam.mViewCells, true);
2892                        break;
2893                case 0:
2894                       
2895                        if (node->IsLeaf())
2896                        {
2897                                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2898                       
2899                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
2900                                {
2901                                        leaf->GetViewCell()->Mail();
2902                                        beam.mViewCells.push_back(leaf->GetViewCell());
2903                                }
2904                        }
2905                        else
2906                        {
2907                                BspInterior *interior = dynamic_cast<BspInterior *>(node);
2908                       
2909                                BspNode *first = interior->GetFront();
2910                                BspNode *second = interior->GetBack();
2911           
2912                                BspNodeGeometry *firstGeom = new BspNodeGeometry();
2913                                BspNodeGeometry *secondGeom = new BspNodeGeometry();
2914
2915                                geom->SplitGeometry(*firstGeom,
2916                                                                        *secondGeom,
2917                                                                        interior->GetPlane(),
2918                                                                        mBox,
2919                                                                        mEpsilon);
2920
2921                                // decide on the order of the nodes
2922                                if (DotProd(beam.mPlanes[0].mNormal,
2923                                        interior->GetPlane().mNormal) > 0)
2924                                {
2925                                        swap(first, second);
2926                                        swap(firstGeom, secondGeom);
2927                                }
2928
2929                                nodeStack.push(bspNodePair(first, firstGeom));
2930                                nodeStack.push(bspNodePair(second, secondGeom));
2931                        }
2932                       
2933                        break;
2934                default:
2935                        // default: cull
2936                        break;
2937                }
2938               
2939                DEL_PTR(geom);
2940               
2941        }
2942
2943        return (int)beam.mViewCells.size();
2944}
2945
2946
2947void VspBspTree::SetViewCellsManager(ViewCellsManager *vcm)
2948{
2949        mViewCellsManager = vcm;
2950}
2951
2952
2953int VspBspTree::CollectMergeCandidates(const vector<BspLeaf *> leaves,
2954                                                                           vector<MergeCandidate> &candidates)
2955{
2956        BspLeaf::NewMail();
2957       
2958        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
2959
2960        int numCandidates = 0;
2961
2962        // find merge candidates and push them into queue
2963        for (it = leaves.begin(); it != it_end; ++ it)
2964        {
2965                BspLeaf *leaf = *it;
2966               
2967                // the same leaves must not be part of two merge candidates
2968                leaf->Mail();
2969               
2970                vector<BspLeaf *> neighbors;
2971                if (0)
2972                        FindNeighbors(leaf, neighbors, true);
2973                else
2974                        FindApproximateNeighbors(leaf, neighbors, true);
2975                vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();
2976
2977                // TODO: test if at least one ray goes from one leaf to the other
2978                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
2979                {
2980                        if ((*nit)->GetViewCell() != leaf->GetViewCell())
2981                        {
2982                                MergeCandidate mc(leaf->GetViewCell(), (*nit)->GetViewCell());
2983                                candidates.push_back(mc);
2984
2985                                ++ numCandidates;
2986                                if ((numCandidates % 1000) == 0)
2987                                {
2988                                        cout << "collected " << numCandidates << " merge candidates" << endl;
2989                                }
2990                        }
2991                }
2992        }
2993
2994        Debug << "merge queue: " << (int)candidates.size() << endl;
2995        Debug << "leaves in queue: " << numCandidates << endl;
2996       
2997
2998        return (int)leaves.size();
2999}
3000
3001
3002int VspBspTree::CollectMergeCandidates(const VssRayContainer &rays,
3003                                                                           vector<MergeCandidate> &candidates)
3004{
3005        ViewCell::NewMail();
3006        long startTime = GetTime();
3007       
3008        map<BspLeaf *, vector<BspLeaf*> > neighborMap;
3009        ViewCellContainer::const_iterator iit;
3010
3011        int numLeaves = 0;
3012       
3013        BspLeaf::NewMail();
3014
3015        for (int i = 0; i < (int)rays.size(); ++ i)
3016        { 
3017                VssRay *ray = rays[i];
3018       
3019                // traverse leaves stored in the rays and compare and
3020                // merge consecutive leaves (i.e., the neighbors in the tree)
3021                if (ray->mViewCells.size() < 2)
3022                        continue;
3023//TODO viewcellhierarchy
3024                iit = ray->mViewCells.begin();
3025                BspViewCell *bspVc = dynamic_cast<BspViewCell *>(*(iit ++));
3026                BspLeaf *leaf = bspVc->mLeaf;
3027               
3028                // traverse intersections
3029                // consecutive leaves are neighbors => add them to queue
3030                for (; iit != ray->mViewCells.end(); ++ iit)
3031                {
3032                        // next pair
3033                        BspLeaf *prevLeaf = leaf;
3034                        bspVc = dynamic_cast<BspViewCell *>(*iit);
3035            leaf = bspVc->mLeaf;
3036
3037                        // view space not valid or same view cell
3038                        if (!leaf->TreeValid() || !prevLeaf->TreeValid() ||
3039                                (leaf->GetViewCell() == prevLeaf->GetViewCell()))
3040                                continue;
3041
3042                vector<BspLeaf *> &neighbors = neighborMap[leaf];
3043                       
3044                        bool found = false;
3045
3046                        // both leaves inserted in queue already =>
3047                        // look if double pair already exists
3048                        if (leaf->Mailed() && prevLeaf->Mailed())
3049                        {
3050                                vector<BspLeaf *>::const_iterator it, it_end = neighbors.end();
3051                               
3052                for (it = neighbors.begin(); !found && (it != it_end); ++ it)
3053                                        if (*it == prevLeaf)
3054                                                found = true; // already in queue
3055                        }
3056               
3057                        if (!found)
3058                        {
3059                                // this pair is not in map yet
3060                                // => insert into the neighbor map and the queue
3061                                neighbors.push_back(prevLeaf);
3062                                neighborMap[prevLeaf].push_back(leaf);
3063
3064                                leaf->Mail();
3065                                prevLeaf->Mail();
3066               
3067                                MergeCandidate mc(leaf->GetViewCell(), prevLeaf->GetViewCell());
3068                               
3069                                candidates.push_back(mc);
3070
3071                                if (((int)candidates.size() % 1000) == 0)
3072                                {
3073                                        cout << "collected " << (int)candidates.size() << " merge candidates" << endl;
3074                                }
3075                        }
3076        }
3077        }
3078
3079        Debug << "neighbormap size: " << (int)neighborMap.size() << endl;
3080        Debug << "merge queue: " << (int)candidates.size() << endl;
3081        Debug << "leaves in queue: " << numLeaves << endl;
3082
3083
3084        //-- collect the leaves which haven't been found by ray casting
3085        if (0)
3086        {
3087                cout << "finding additional merge candidates using geometry" << endl;
3088                vector<BspLeaf *> leaves;
3089                CollectLeaves(leaves, true);
3090                Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
3091                CollectMergeCandidates(leaves, candidates);
3092        }
3093
3094        return numLeaves;
3095}
3096
3097
3098
3099
3100ViewCell *VspBspTree::GetViewCell(const Vector3 &point)
3101{
3102  if (mRoot == NULL)
3103        return NULL;
3104 
3105  stack<BspNode *> nodeStack;
3106  nodeStack.push(mRoot);
3107 
3108  ViewCell *viewcell = NULL;
3109 
3110  while (!nodeStack.empty())  {
3111        BspNode *node = nodeStack.top();
3112        nodeStack.pop();
3113       
3114        if (node->IsLeaf()) {
3115          viewcell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
3116          break;
3117        } else {
3118         
3119          BspInterior *interior = dynamic_cast<BspInterior *>(node);
3120               
3121          // random decision
3122          if (interior->GetPlane().Side(point) < 0)
3123                nodeStack.push(interior->GetBack());
3124          else
3125                nodeStack.push(interior->GetFront());
3126        }
3127  }
3128 
3129  return viewcell;
3130}
3131
3132
3133bool VspBspTree::ViewPointValid(const Vector3 &viewPoint) const
3134{
3135        BspNode *node = mRoot;
3136
3137        while (1)
3138        {
3139                // early exit
3140                if (node->TreeValid())
3141                        return true;
3142
3143                if (node->IsLeaf())
3144                        return false;
3145                       
3146                BspInterior *in = dynamic_cast<BspInterior *>(node);
3147                                       
3148                if (in->GetPlane().Side(viewPoint) <= 0)
3149                {
3150                        node = in->GetBack();
3151                }
3152                else
3153                {
3154                        node = in->GetFront();
3155                }
3156        }
3157
3158        // should never come here
3159        return false;
3160}
3161
3162
3163void VspBspTree::PropagateUpValidity(BspNode *node)
3164{
3165        const bool isValid = node->TreeValid();
3166
3167        // propagative up invalid flag until only invalid nodes exist over this node
3168        if (!isValid)
3169        {
3170                while (!node->IsRoot() && node->GetParent()->TreeValid())
3171                {
3172                        node = node->GetParent();
3173                        node->SetTreeValid(false);
3174                }
3175        }
3176        else
3177        {
3178                // propagative up valid flag until one of the subtrees is invalid
3179                while (!node->IsRoot() && !node->TreeValid())
3180                {
3181            node = node->GetParent();
3182                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
3183                       
3184                        // the parent is valid iff both leaves are valid
3185                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
3186                                                           interior->GetFront()->TreeValid());
3187                }
3188        }
3189}
3190
3191
3192bool VspBspTree::Export(ofstream &stream)
3193{
3194        ExportNode(mRoot, stream);
3195
3196        return true;
3197}
3198
3199
3200void VspBspTree::ExportNode(BspNode *node, ofstream &stream)
3201{
3202        if (node->IsLeaf())
3203        {
3204                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
3205                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
3206
3207                int id = -1;
3208                if (viewCell != mOutOfBoundsCell)
3209                        id = viewCell->GetId();
3210
3211                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
3212        }
3213        else
3214        {
3215                BspInterior *interior = dynamic_cast<BspInterior *>(node);
3216       
3217                Plane3 plane = interior->GetPlane();
3218                stream << "<Interior plane=\"" << plane.mNormal.x << " "
3219                           << plane.mNormal.y << " " << plane.mNormal.z << " "
3220                           << plane.mD << "\">" << endl;
3221
3222                ExportNode(interior->GetBack(), stream);
3223                ExportNode(interior->GetFront(), stream);
3224
3225                stream << "</Interior>" << endl;
3226        }
3227}
Note: See TracBrowser for help on using the repository browser.