source: GTP/trunk/Lib/Vis/Preprocessing/src/VspOspTree.cpp @ 1002

Revision 1002, 94.5 KB checked in by mattausch, 18 years ago (diff)

debug run: fixing memory holes

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