source: GTP/trunk/Lib/Vis/Preprocessing/src/FromPointVisibilityTree.cpp @ 1715

Revision 1715, 97.3 KB checked in by bittner, 18 years ago (diff)

new visibility filter support

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