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

Revision 882, 93.9 KB checked in by mattausch, 18 years ago (diff)

added new active view cell concept

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