source: trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp @ 607

Revision 607, 72.5 KB checked in by mattausch, 18 years ago (diff)

added method for associating spatial hierarchy leaf with view cell

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