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

Revision 582, 69.7 KB checked in by mattausch, 18 years ago (diff)

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