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

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