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

Revision 445, 46.3 KB checked in by mattausch, 19 years ago (diff)

Added VspBspTree? functionality:
This is a view space partitioning specialised BSPtree.
There are no polygon splits, but we split the sample rays.
The candidates for the next split plane are evaluated only
by checking the sampled visibility information.
The polygons are employed merely as candidates for the next split planes.
Line 
1#include "Plane3.h"
2#include "VspBspTree.h"
3#include "Mesh.h"
4#include "common.h"
5#include "ViewCell.h"
6#include "Environment.h"
7#include "Polygon3.h"
8#include "Ray.h"
9#include "AxisAlignedBox3.h"
10#include <stack>
11#include <time.h>
12#include <iomanip>
13#include "Exporter.h"
14#include "Plane3.h"
15
16//-- static members
17/** Evaluates split plane classification with respect to the plane's
18        contribution for a minimum number of ray splits.
19*/
20const float VspBspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
21/** Evaluates split plane classification with respect to the plane's
22        contribution for balanced rays.
23*/
24const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
25
26int VspBspNode::sMailId = 1;
27
28int VspBspTree::sFrontId = 0;
29int VspBspTree::sBackId = 0;
30int VspBspTree::sFrontAndBackId = 0;
31
32/****************************************************************/
33/*                  class VspBspNode implementation                */
34/****************************************************************/
35
36VspBspNode::VspBspNode():
37mParent(NULL)
38{}
39
40VspBspNode::VspBspNode(VspBspInterior *parent):
41mParent(parent)
42{}
43
44
45bool VspBspNode::IsRoot() const
46{
47        return mParent == NULL;
48}
49
50VspBspInterior *VspBspNode::GetParent()
51{
52        return mParent;
53}
54
55void VspBspNode::SetParent(VspBspInterior *parent)
56{
57        mParent = parent;
58}
59
60/****************************************************************/
61/*              class VspBspInterior implementation                */
62/****************************************************************/
63
64
65VspBspInterior::VspBspInterior(const Plane3 &plane):
66mPlane(plane), mFront(NULL), mBack(NULL)
67{}
68
69VspBspInterior::~VspBspInterior()
70{
71        DEL_PTR(mFront);
72        DEL_PTR(mBack);
73}
74
75bool VspBspInterior::IsLeaf() const
76{
77        return false;
78}
79
80VspBspNode *VspBspInterior::GetBack()
81{
82        return mBack;
83}
84
85VspBspNode *VspBspInterior::GetFront()
86{
87        return mFront;
88}
89
90Plane3 *VspBspInterior::GetPlane()
91{
92        return &mPlane;
93}
94
95void VspBspInterior::ReplaceChildLink(VspBspNode *oldChild, VspBspNode *newChild)
96{
97        if (mBack == oldChild)
98                mBack = newChild;
99        else
100                mFront = newChild;
101}
102
103void VspBspInterior::SetupChildLinks(VspBspNode *b, VspBspNode *f)
104{
105    mBack = b;
106    mFront = f;
107}
108
109int VspBspInterior::SplitPolygons(PolygonContainer &polys,
110                                                                  PolygonContainer &frontPolys,
111                                                                  PolygonContainer &backPolys,
112                                                                  PolygonContainer &coincident)
113{
114        int splits = 0;
115
116        while (!polys.empty())
117        {
118                Polygon3 *poly = polys.back();
119                polys.pop_back();
120
121                // classify polygon
122                const int cf = poly->ClassifyPlane(mPlane);
123
124                switch (cf)
125                {
126                        case Polygon3::COINCIDENT:
127                                coincident.push_back(poly);
128                                break;                 
129                        case Polygon3::FRONT_SIDE:     
130                                frontPolys.push_back(poly);
131                                break;
132                        case Polygon3::BACK_SIDE:
133                                backPolys.push_back(poly);
134                                break;
135                        case Polygon3::SPLIT:
136                                backPolys.push_back(poly);
137                                frontPolys.push_back(poly);
138                                ++ splits;
139                                break;
140                        default:
141                Debug << "SHOULD NEVER COME HERE\n";
142                                break;
143                }
144        }
145
146        return splits;
147}
148
149/****************************************************************/
150/*                  class VspBspLeaf implementation                */
151/****************************************************************/
152VspBspLeaf::VspBspLeaf(): mViewCell(NULL)
153{
154}
155
156VspBspLeaf::VspBspLeaf(BspViewCell *viewCell):
157mViewCell(viewCell)
158{
159}
160
161VspBspLeaf::VspBspLeaf(VspBspInterior *parent):
162VspBspNode(parent), mViewCell(NULL)
163{}
164
165VspBspLeaf::VspBspLeaf(VspBspInterior *parent, BspViewCell *viewCell):
166VspBspNode(parent), mViewCell(viewCell)
167{
168}
169
170BspViewCell *VspBspLeaf::GetViewCell() const
171{
172        return mViewCell;
173}
174
175void VspBspLeaf::SetViewCell(BspViewCell *viewCell)
176{
177        mViewCell = viewCell;
178}
179
180bool VspBspLeaf::IsLeaf() const
181{
182        return true;
183}
184
185void VspBspLeaf::AddToPvs(const RayInfoContainer &rays,
186                                                  int &sampleContributions,
187                                                  int &contributingSamples,
188                                                  bool storeLeavesWithRays)
189{
190        sampleContributions = 0;
191        contributingSamples = 0;
192
193    RayInfoContainer::const_iterator it, it_end = rays.end();
194
195        // add contributions from samples to the PVS
196        for (it = rays.begin(); it != it_end; ++ it)
197        {
198                int contribution = 0;
199                VssRay *ray = (*it).mRay;
200                       
201                if (ray->mTerminationObject)
202                        contribution += mViewCell->GetPvs().AddSample(ray->mTerminationObject);
203               
204                if (ray->mOriginObject)
205                        contribution += mViewCell->GetPvs().AddSample(ray->mOriginObject);
206
207                if (contribution > 0)
208                {
209                        sampleContributions += contribution;
210                        ++ contributingSamples;
211                }
212        }
213}
214
215
216/****************************************************************/
217/*                  class VspBspTree implementation             */
218/****************************************************************/
219
220VspBspTree::VspBspTree():
221mRoot(NULL),
222mPvsUseArea(true)
223{
224        mRootCell = new BspViewCell();
225
226        Randomize(); // initialise random generator for heuristics
227
228        //-- termination criteria for autopartition
229        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
230        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
231        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
232        environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);     
233        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
234        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
235
236        //-- factors for bsp tree split plane heuristics
237        environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
238        environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
239        environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
240
241        //-- termination criteria for axis aligned split
242        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAaCtDivCi);
243        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
244       
245        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
246                                                         mTermMinRaysForAxisAligned);
247       
248        //-- partition criteria
249        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
250        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
251        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
252        environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mSplitBorder);
253
254        environment->GetFloatValue("VspBspTree.Construction.sideTolerance", Vector3::sDistTolerance);
255        Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
256
257    Debug << "BSP max depth: " << mTermMaxDepth << endl;
258        Debug << "BSP min PVS: " << mTermMinPvs << endl;
259        Debug << "BSP min area: " << mTermMinArea << endl;
260        Debug << "BSP min rays: " << mTermMinRays << endl;
261        Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
262        Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;
263
264        Debug << "Split plane strategy: ";
265        if (mSplitPlaneStrategy & RANDOM_POLYGON)
266                Debug << "random polygon ";
267        if (mSplitPlaneStrategy & AXIS_ALIGNED)
268                Debug << "axis aligned ";
269        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
270                Debug << "least ray splits ";
271        if (mSplitPlaneStrategy & BALANCED_RAYS)
272                Debug << "balanced rays ";
273        if (mSplitPlaneStrategy & PVS)
274                Debug << "pvs";
275
276        Debug << endl;
277}
278
279
280const VspBspTreeStatistics &VspBspTree::GetStatistics() const
281{
282        return mStat;
283}
284
285void VspBspTreeStatistics::Print(ostream &app) const
286{
287        app << "===== VspBspTree statistics ===============\n";
288
289        app << setprecision(4);
290
291        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
292
293        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
294
295        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
296
297        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
298
299        app << "#N_SPLITS ( Number of splits )\n" << splits << "\n";
300
301        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
302                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
303
304        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maximum depth )\n"
305                <<      maxDepthNodes * 100 / (double)Leaves() << endl;
306
307        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n"
308                << minPvsNodes * 100 / (double)Leaves() << endl;
309
310        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n"
311                <<      minRaysNodes * 100 / (double)Leaves() << endl;
312
313        app << "#N_PMINAREALEAVES  ( Percentage of leaves with mininum area )\n"
314                << minAreaNodes * 100 / (double)Leaves() << endl;
315
316        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"
317                <<      maxRayContribNodes * 100 / (double)Leaves() << endl;
318
319        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
320
321        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
322
323        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
324
325        app << "#N_INPUT_POLYGONS (number of input polygons )\n" <<     polys << endl;
326
327        //app << "#N_PVS: " << pvs << endl;
328
329        app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
330                 (polys + splits) / (double)polys << endl;
331       
332        app << "===== END OF VspBspTree statistics ==========\n";
333}
334
335
336VspBspTree::~VspBspTree()
337{
338        DEL_PTR(mRoot);
339        DEL_PTR(mRootCell);
340}
341
342int VspBspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
343{
344        FaceContainer::const_iterator fi;
345       
346        // copy the face data to polygons
347        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
348        {
349                Polygon3 *poly = new Polygon3((*fi), mesh);
350               
351                if (poly->Valid())
352                {
353                        poly->mParent = parent; // set parent intersectable
354                        polys.push_back(poly);
355                }
356                else
357                        DEL_PTR(poly);
358        }
359        return (int)mesh->mFaces.size();
360}
361
362int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells,
363                                                          PolygonContainer &polys,
364                                                          int maxObjects)
365{
366        int limit = (maxObjects > 0) ?
367                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
368 
369        int polysSize = 0;
370
371        for (int i = 0; i < limit; ++ i)
372        {
373                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
374                {
375                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
376                        polysSize += AddMeshToPolygons(viewCells[i]->GetMesh(), polys, viewCells[i]);
377                }
378        }
379
380        return polysSize;
381}
382
383int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects, PolygonContainer &polys, int maxObjects)
384{
385        int limit = (maxObjects > 0) ? Min((int)objects.size(), maxObjects) : (int)objects.size();
386 
387        for (int i = 0; i < limit; ++i)
388        {
389                Intersectable *object = objects[i];//*it;
390                Mesh *mesh = NULL;
391
392                switch (object->Type()) // extract the meshes
393                {
394                case Intersectable::MESH_INSTANCE:
395                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
396                        break;
397                case Intersectable::VIEW_CELL:
398                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
399                        break;
400                        // TODO: handle transformed mesh instances
401                default:
402                        Debug << "intersectable type not supported" << endl;
403                        break;
404                }
405               
406        if (mesh) // copy the mesh data to polygons
407                {
408                        mBox.Include(object->GetBox()); // add to BSP tree aabb
409                        AddMeshToPolygons(mesh, polys, mRootCell);
410                }
411        }
412
413        return (int)polys.size();
414}
415
416void VspBspTree::Construct(const VssRayContainer &sampleRays)
417{
418    mStat.nodes = 1;
419        mBox.Initialize();      // initialise BSP tree bounding box
420       
421        PolygonContainer polys;
422        RayInfoContainer *rays = new RayInfoContainer();
423
424        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
425
426        long startTime = GetTime();
427
428        Debug << "**** Extracting polygons from rays ****\n";
429
430        //-- extract polygons intersected by the rays
431        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
432        {
433                VssRay *ray = *rit;
434       
435                if (ray->mTerminationObject)
436                {
437                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
438                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
439                }
440
441                if (ray->mOriginObject)
442                {
443                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
444                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
445                }
446        }
447
448        //-- store rays
449        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
450        {
451                rays->push_back(RayInfo(*rit));
452        }
453
454        mStat.polys = (int)polys.size();
455
456        Debug << "**** Finished polygon extraction ****" << endl;
457        Debug << (int)polys.size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
458        Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
459
460        Construct(new PolygonContainer(polys), rays);
461
462        // clean up polygons because they are not deleted in the tree
463        CLEAR_CONTAINER(polys);
464}
465
466void VspBspTree::Construct(PolygonContainer *polys, RayInfoContainer *rays)
467{
468        std::stack<VspBspTraversalData> tStack;
469
470        mRoot = new VspBspLeaf();
471
472        VspBspNodeGeometry *cell = new VspBspNodeGeometry();
473        ConstructGeometry(mRoot, *cell);
474
475        VspBspTraversalData tData(mRoot, polys, 0, mRootCell, rays,
476                                                   ComputePvsSize(*rays), cell->GetArea(), cell);
477
478        tStack.push(tData);
479
480        mStat.Start();
481        cout << "Contructing bsp tree ... ";
482        long startTime = GetTime();
483        while (!tStack.empty())
484        {
485                tData = tStack.top();
486
487            tStack.pop();
488
489                // subdivide leaf node
490                VspBspNode *r = Subdivide(tStack, tData);
491
492                if (r == mRoot)
493                        Debug << "BSP tree construction time spent at root: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
494        }
495
496        cout << "finished\n";
497
498        mStat.Stop();
499}
500
501bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data) const
502{
503        return
504                (((int)data.mRays->size() <= mTermMinRays) ||
505                 (data.mPvs <= mTermMinPvs) ||
506                 (data.mArea <= mTermMinArea) ||
507                 (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
508                 (data.mDepth >= mTermMaxDepth));
509}
510
511VspBspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack, VspBspTraversalData &tData)
512{
513        //-- terminate traversal 
514        if (TerminationCriteriaMet(tData))             
515        {
516                VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(tData.mNode);
517       
518                BspViewCell *viewCell = new BspViewCell();
519               
520                leaf->SetViewCell(viewCell);
521                //TODO viewCell->mVspBspLeaves.push_back(leaf);
522
523                //-- add pvs
524                if (viewCell != mRootCell)
525                {
526                        int conSamp = 0, sampCon = 0;
527                        leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
528                       
529                        mStat.contributingSamples += conSamp;
530                        mStat.sampleContributions += sampCon;
531                }
532
533                EvaluateLeafStats(tData);
534               
535                //-- clean up
536
537                DEL_PTR(tData.mPolygons);
538                DEL_PTR(tData.mRays);
539                DEL_PTR(tData.mGeometry);
540               
541                return leaf;
542        }
543
544        //-- continue subdivision
545        PolygonContainer coincident;
546       
547        VspBspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell,
548                                                                new RayInfoContainer(), 0, 0, new VspBspNodeGeometry());
549        VspBspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell,
550                                                           new RayInfoContainer(), 0, 0, new VspBspNodeGeometry());
551
552        // create new interior node and two leaf nodes
553        VspBspInterior *interior =
554                SubdivideNode(tData, tFrontData, tBackData, coincident);
555
556#ifdef _DEBUG   
557        if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
558        {       for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
559                        Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
560                Debug << endl;}
561#endif
562
563        // don't need coincident polygons anymory
564        CLEAR_CONTAINER(coincident);
565
566        // push the children on the stack
567        tStack.push(tFrontData);
568        tStack.push(tBackData);
569
570        // cleanup
571        DEL_PTR(tData.mNode);
572        DEL_PTR(tData.mPolygons);
573        DEL_PTR(tData.mRays);
574        DEL_PTR(tData.mGeometry);               
575
576        return interior;
577}
578
579VspBspInterior *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
580                                                                        VspBspTraversalData &frontData,
581                                                                        VspBspTraversalData &backData,
582                                                                        PolygonContainer &coincident)
583{
584        mStat.nodes += 2;
585       
586        VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(tData.mNode);
587        // select subdivision plane
588        VspBspInterior *interior =
589                new VspBspInterior(SelectPlane(leaf, tData));
590
591#ifdef _DEBUG
592        Debug << interior << endl;
593#endif
594       
595        // subdivide rays into front and back rays
596        SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);
597       
598        // subdivide polygons with plane
599        mStat.splits += interior->SplitPolygons(*tData.mPolygons,
600                                                    *frontData.mPolygons,
601                                                                                        *backData.mPolygons,
602                                                                                        coincident);
603
604    // compute pvs
605        frontData.mPvs = ComputePvsSize(*frontData.mRays);
606        backData.mPvs = ComputePvsSize(*backData.mRays);
607
608        // split geometry and compute area
609        if (1)
610        {
611                tData.mGeometry->SplitGeometry(*frontData.mGeometry,
612                                                                           *backData.mGeometry,
613                                                                           *this,
614                                                                           interior->mPlane);
615       
616               
617                frontData.mArea = frontData.mGeometry->GetArea();
618                backData.mArea = backData.mGeometry->GetArea();
619        }
620
621        // compute accumulated ray length
622        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
623        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
624
625        //-- create front and back leaf
626
627        VspBspInterior *parent = leaf->GetParent();
628
629        // replace a link from node's parent
630        if (!leaf->IsRoot())
631        {
632                parent->ReplaceChildLink(leaf, interior);
633                interior->SetParent(parent);
634        }
635        else // new root
636        {
637                mRoot = interior;
638        }
639
640        // and setup child links
641        interior->SetupChildLinks(new VspBspLeaf(interior), new VspBspLeaf(interior));
642       
643        frontData.mNode = interior->mFront;
644        backData.mNode = interior->mBack;
645       
646        //DEL_PTR(leaf);
647        return interior;
648}
649
650void VspBspTree::SortSplitCandidates(const PolygonContainer &polys,
651                                                                  const int axis,
652                                                                  vector<SortableEntry> &splitCandidates) const
653{
654  splitCandidates.clear();
655 
656  int requestedSize = 2 * (int)polys.size();
657  // creates a sorted split candidates array 
658  splitCandidates.reserve(requestedSize);
659 
660  PolygonContainer::const_iterator it, it_end = polys.end();
661
662  AxisAlignedBox3 box;
663
664  // insert all queries
665  for(it = polys.begin(); it != it_end; ++ it)
666  {
667          box.Initialize();
668          box.Include(*(*it));
669         
670          splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
671      splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
672  }
673 
674  stable_sort(splitCandidates.begin(), splitCandidates.end());
675}
676
677
678float VspBspTree::BestCostRatio(const PolygonContainer &polys,
679                                                         const AxisAlignedBox3 &box,
680                                                         const int axis,
681                                                         float &position,
682                                                         int &objectsBack,
683                                                         int &objectsFront) const
684{
685        vector<SortableEntry> splitCandidates;
686
687        SortSplitCandidates(polys, axis, splitCandidates);
688       
689        // go through the lists, count the number of objects left and right
690        // and evaluate the following cost funcion:
691        // C = ct_div_ci  + (ol + or)/queries
692       
693        int objectsLeft = 0, objectsRight = (int)polys.size();
694       
695        float minBox = box.Min(axis);
696        float maxBox = box.Max(axis);
697        float boxArea = box.SurfaceArea();
698 
699        float minBand = minBox + mSplitBorder * (maxBox - minBox);
700        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
701       
702        float minSum = 1e20f;
703        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
704
705        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci)
706        {
707                switch ((*ci).type)
708                {
709                        case SortableEntry::POLY_MIN:
710                                ++ objectsLeft;
711                                break;
712                        case SortableEntry::POLY_MAX:
713                            -- objectsRight;
714                                break;
715                        default:
716                                break;
717                }
718               
719                if ((*ci).value > minBand && (*ci).value < maxBand)
720                {
721                        AxisAlignedBox3 lbox = box;
722                        AxisAlignedBox3 rbox = box;
723                        lbox.SetMax(axis, (*ci).value);
724                        rbox.SetMin(axis, (*ci).value);
725
726                        float sum = objectsLeft * lbox.SurfaceArea() +
727                                                objectsRight * rbox.SurfaceArea();
728     
729                        if (sum < minSum)
730                        {
731                                minSum = sum;
732                                position = (*ci).value;
733
734                                objectsBack = objectsLeft;
735                                objectsFront = objectsRight;
736                        }
737                }
738        }
739 
740        float oldCost = (float)polys.size();
741        float newCost = mAaCtDivCi + minSum / boxArea;
742        float ratio = newCost / oldCost;
743
744
745#if 0
746  Debug << "====================" << endl;
747  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
748      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
749#endif
750  return ratio;
751}
752
753bool VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
754                                                                         const PolygonContainer &polys) const
755{
756        AxisAlignedBox3 box;
757        box.Initialize();
758       
759        // create bounding box of region
760        Polygon3::IncludeInBox(polys, box);
761       
762        int objectsBack = 0, objectsFront = 0;
763        int axis = 0;
764        float costRatio = MAX_FLOAT;
765        Vector3 position;
766
767        //-- area subdivision
768        for (int i = 0; i < 3; ++ i)
769        {
770                float p = 0;
771                float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
772               
773                if (r < costRatio)
774                {
775                        costRatio = r;
776                        axis = i;
777                        position = p;
778                }
779        }
780       
781        if (costRatio >= mMaxCostRatio)
782                return false;
783
784        Vector3 norm(0,0,0); norm[axis] = 1.0f;
785        plane = Plane3(norm, position);
786
787        return true;
788}
789
790Plane3 VspBspTree::SelectPlane(VspBspLeaf *leaf, VspBspTraversalData &data)
791{
792        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
793                ((int)data.mRays->size() > mTermMinRaysForAxisAligned))
794        {
795                Plane3 plane;
796                if (SelectAxisAlignedPlane(plane, *data.mPolygons)) // TODO: should be rays
797                        return plane;
798        }
799
800        // simplest strategy: just take next polygon
801        if (mSplitPlaneStrategy & RANDOM_POLYGON)
802        {
803        if (!data.mPolygons->empty())
804                {
805                        Polygon3 *nextPoly = (*data.mPolygons)[Random((int)data.mPolygons->size())];
806                        return nextPoly->GetSupportingPlane();
807                }
808                else
809                {
810                        //-- choose plane on midpoint of a ray
811                        const int candidateIdx = Random((int)data.mRays->size());
812                                                               
813                        const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
814                        const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();
815
816                        const Vector3 pt = (maxPt + minPt) * 0.5;
817
818                        const Vector3 normal =  (*data.mRays)[candidateIdx].mRay->GetDir();
819                       
820                        return Plane3(normal, pt);
821                }
822
823                return Plane3();
824        }
825
826        // use heuristics to find appropriate plane
827        return SelectPlaneHeuristics(leaf, data);
828}
829
830Plane3 VspBspTree::SelectPlaneHeuristics(VspBspLeaf *leaf, VspBspTraversalData &data)
831{
832        float lowestCost = MAX_FLOAT;
833        Plane3 bestPlane;
834        Plane3 plane;
835
836        int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
837       
838        int candidateIdx = limit;
839
840        for (int i = 0; i < limit; ++ i)
841        {
842                candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
843               
844                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
845                // evaluate current candidate
846                const float candidateCost =
847                        SplitPlaneCost(poly->GetSupportingPlane(), data);
848
849                if (candidateCost < lowestCost)
850                {
851                        bestPlane = poly->GetSupportingPlane();
852                        lowestCost = candidateCost;
853                }
854        }
855       
856        //Debug << "lowest: " << lowestCost << endl;
857
858        //-- choose candidate planes extracted from rays
859        // we currently use different two methods chosen with
860        // equal probability
861       
862        RayInfoContainer *rays = data.mRays;
863        // take 3 ray endpoints, where two are minimum and one a maximum
864        // point or the other way round
865        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
866        {
867                candidateIdx = Random((int)rays->size());
868       
869                RayInfo rayInf = (*rays)[candidateIdx];
870                const Vector3 minPt = rayInf.ExtrapOrigin();
871                const Vector3 maxPt = rayInf.ExtrapTermination();
872
873                const Vector3 pt = (maxPt + minPt) * 0.5;
874                const Vector3 normal = Normalize(rayInf.mRay->GetDir());
875
876                plane = Plane3(normal, pt);
877       
878                const float candidateCost = SplitPlaneCost(plane, data);
879
880                if (candidateCost < lowestCost)
881                {
882                        bestPlane = plane;     
883                        lowestCost = candidateCost;
884                }
885        }
886
887        // take plane normal as plane normal and the midpoint of the ray.
888        // PROBLEM: does not resemble any point where visibility is likely to change
889        //Debug << "lowest: " << lowestCost << endl;
890        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
891        {
892                Vector3 pt[3];
893                int idx[3];
894                int cmaxT = 0;
895                int cminT = 0;
896                bool chooseMin = false;
897
898                for (int j = 0; j < 3; j ++)
899                {
900                        idx[j] = Random((int)rays->size() * 2);
901                               
902                        if (idx[j] >= (int)rays->size())
903                        {
904                                idx[j] -= (int)rays->size();
905                               
906                                chooseMin = (cminT < 2);
907                        }
908                        else
909                                chooseMin = (cmaxT < 2);
910
911                        RayInfo rayInf = (*rays)[idx[j]];
912                        pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
913                }       
914                       
915                plane = Plane3(pt[0], pt[1], pt[2]);
916
917                const float candidateCost = SplitPlaneCost(plane, data);
918
919                if (candidateCost < lowestCost)
920                {
921                        //Debug << "choose ray plane 2: " << candidateCost << endl;
922                        bestPlane = plane;
923                       
924                        lowestCost = candidateCost;
925                }
926        }       
927
928#ifdef _DEBUG
929        Debug << "plane lowest cost: " << lowestCost << endl;
930#endif
931        return bestPlane;
932}
933
934int VspBspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
935{
936        const int candidateIdx = Random(currentIdx --);
937
938        // swap candidates to avoid testing same plane 2 times
939        std::swap(polys[currentIdx], polys[candidateIdx]);
940       
941        return currentIdx;
942        //return Random((int)polys.size());
943}
944
945
946bool VspBspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
947{
948        maxT = 1e6;
949        minT = 0;
950       
951        // test with tree bounding box
952        if (!mBox.GetMinMaxT(ray, &minT, &maxT))
953                return false;
954
955        if (minT < 0) // start ray from origin
956                minT = 0;
957
958        // bound ray or line segment
959        if ((ray.GetType() == Ray::LOCAL_RAY) &&
960            !ray.intersections.empty() &&
961                (ray.intersections[0].mT <= maxT))
962        {
963                maxT = ray.intersections[0].mT;
964        }
965
966        return true;
967}
968
969inline void VspBspTree::GenerateUniqueIdsForPvs()
970{
971        Intersectable::NewMail(); sBackId = ViewCell::sMailId;
972        Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
973        Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
974}
975
976float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane,
977                                                                 const VspBspTraversalData &data)
978{
979        float val = 0;
980
981        float sumBalancedRays = 0;
982        float sumRaySplits = 0;
983       
984        int backId = 0;
985        int frontId = 0;
986        int frontAndBackId = 0;
987
988        int frontPvs = 0;
989        int backPvs = 0;
990
991        // probability that view point lies in child
992        float pOverall = 0;
993        float pFront = 0;
994        float pBack = 0;
995
996        const bool pvsUseLen = false;
997
998        if (mSplitPlaneStrategy & PVS)
999        {
1000                // create unique ids for pvs heuristics
1001                GenerateUniqueIdsForPvs();
1002       
1003                if (mPvsUseArea) // use front and back cell areas to approximate volume
1004                {       
1005                        // construct child geometry with regard to the candidate split plane
1006                        VspBspNodeGeometry frontCell;
1007                        VspBspNodeGeometry backCell;
1008
1009                        data.mGeometry->SplitGeometry(frontCell, backCell, *this, candidatePlane);
1010               
1011                        pFront = frontCell.GetArea();
1012                        pBack = backCell.GetArea();
1013
1014                        pOverall = data.mArea;
1015                }
1016        }
1017                       
1018        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1019
1020        for (rit = data.mRays->begin(); rit != data.mRays->end(); ++ rit)
1021        {
1022                VssRay *ray = (*rit).mRay;
1023                const int cf = (*rit).ComputeRayIntersection(candidatePlane, (*rit).mRay->mT);
1024
1025                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1026                {
1027                        sumBalancedRays += cf;
1028                }
1029               
1030                if (mSplitPlaneStrategy & BALANCED_RAYS)
1031                {
1032                        if (cf == 0)
1033                                ++ sumRaySplits;
1034                }
1035
1036                if (mSplitPlaneStrategy & PVS)
1037                {
1038                        // in case the ray intersects an object
1039                        // assure that we only count the object
1040                        // once for the front and once for the back side of the plane
1041                       
1042                        // add the termination object
1043                        AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
1044                       
1045                        // add the source object
1046                        AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
1047                       
1048                        // use number or length of rays to approximate volume
1049                        if (!mPvsUseArea)
1050                        {
1051                                float len = 1;
1052
1053                                if (pvsUseLen) // use length of rays
1054                                        len = (*rit).SqrSegmentLength();
1055                       
1056                                pOverall += len;
1057
1058                                if (cf == 1)
1059                                        pFront += len;
1060                                if (cf == -1)
1061                                        pBack += len;
1062                                if (cf == 0)
1063                                {
1064                                        // use length of rays to approximate volume
1065                                        if (pvsUseLen)
1066                                        {
1067                                                float newLen = len *
1068                                                        ((*rit).GetMaxT() - (*rit).mRay->mT) /
1069                                                        ((*rit).GetMaxT() - (*rit).GetMinT());
1070               
1071                                                if (candidatePlane.Side((*rit).ExtrapOrigin()) <= 0)
1072                                                {
1073                                                        pBack += newLen;
1074                                                        pFront += len - newLen;
1075                                                }
1076                                                else
1077                                                {
1078                                                        pFront += newLen;
1079                                                        pBack += len - newLen;
1080                                                }
1081                                        }
1082                                        else
1083                                        {
1084                                                ++ pFront;
1085                                                ++ pBack;
1086                                        }
1087                                }
1088                        }
1089                }
1090        }
1091
1092        const float raysSize = (float)data.mRays->size() + Limits::Small;
1093
1094        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1095                val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
1096
1097        if (mSplitPlaneStrategy & BALANCED_RAYS)
1098                val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
1099
1100        const float denom = pOverall * (float)data.mPvs * 2.0f + Limits::Small;
1101
1102        if (mSplitPlaneStrategy & PVS)
1103        {
1104                val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
1105
1106                // give penalty to unbalanced split
1107                if (0)
1108                if (((pFront * 0.2 + Limits::Small) > pBack) ||
1109                        (pFront < (pBack * 0.2 + Limits::Small)))
1110                        val += 0.5;
1111        }
1112
1113#ifdef _DEBUG
1114        Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
1115                  << " frontpvs: " << frontPvs << " pFront: " << pFront
1116                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
1117#endif
1118        return val;
1119}
1120
1121void VspBspTree::AddObjToPvs(Intersectable *obj,
1122                                                  const int cf,
1123                                                  int &frontPvs,
1124                                                  int &backPvs) const
1125{
1126        if (!obj)
1127                return;
1128        // TODO: does this really belong to no pvs?
1129        //if (cf == Ray::COINCIDENT) return;
1130
1131        // object belongs to both PVS
1132        if (cf >= 0)
1133        {
1134                if ((obj->mMailbox != sFrontId) &&
1135                        (obj->mMailbox != sFrontAndBackId))
1136                {
1137                        ++ frontPvs;
1138
1139                        if (obj->mMailbox == sBackId)
1140                                obj->mMailbox = sFrontAndBackId;       
1141                        else
1142                                obj->mMailbox = sFrontId;                                                               
1143                }
1144        }
1145       
1146        if (cf <= 0)
1147        {
1148                if ((obj->mMailbox != sBackId) &&
1149                        (obj->mMailbox != sFrontAndBackId))
1150                {
1151                        ++ backPvs;
1152
1153                        if (obj->mMailbox == sFrontId)
1154                                obj->mMailbox = sFrontAndBackId;
1155                        else
1156                                obj->mMailbox = sBackId;                               
1157                }
1158        }
1159}
1160
1161void VspBspTree::CollectLeaves(vector<VspBspLeaf *> &leaves) const
1162{
1163        stack<VspBspNode *> nodeStack;
1164        nodeStack.push(mRoot);
1165 
1166        while (!nodeStack.empty())
1167        {
1168                VspBspNode *node = nodeStack.top();
1169   
1170                nodeStack.pop();
1171   
1172                if (node->IsLeaf())
1173                {
1174                        VspBspLeaf *leaf = (VspBspLeaf *)node;         
1175                        leaves.push_back(leaf);
1176                }
1177                else
1178                {
1179                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
1180
1181                        nodeStack.push(interior->GetBack());
1182                        nodeStack.push(interior->GetFront());
1183                }
1184        }
1185}
1186
1187AxisAlignedBox3 VspBspTree::GetBoundingBox() const
1188{
1189        return mBox;
1190}
1191
1192VspBspNode *VspBspTree::GetRoot() const
1193{
1194        return mRoot;
1195}
1196
1197void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
1198{
1199        // the node became a leaf -> evaluate stats for leafs
1200        VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(data.mNode);
1201
1202        // store maximal and minimal depth
1203        if (data.mDepth > mStat.maxDepth)
1204                mStat.maxDepth = data.mDepth;
1205
1206        if (data.mDepth < mStat.minDepth)
1207                mStat.minDepth = data.mDepth;
1208
1209        // accumulate depth to compute average depth
1210        mStat.accumDepth += data.mDepth;
1211       
1212
1213        if (data.mDepth >= mTermMaxDepth)
1214                ++ mStat.maxDepthNodes;
1215
1216        if (data.mPvs < mTermMinPvs)
1217                ++ mStat.minPvsNodes;
1218
1219        if ((int)data.mRays->size() < mTermMinRays)
1220                ++ mStat.minRaysNodes;
1221
1222        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1223                ++ mStat.maxRayContribNodes;
1224       
1225        if (data.mGeometry->GetArea() <= mTermMinArea)
1226                ++ mStat.minAreaNodes;
1227
1228#ifdef _DEBUG
1229        Debug << "BSP stats: "
1230                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1231                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1232                  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
1233                  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << mTermMinPolys << "), "
1234                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1235                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
1236                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1237#endif
1238}
1239
1240int VspBspTree::CastRay(Ray &ray)
1241{
1242        int hits = 0;
1243 
1244        stack<VspBspRayTraversalData> tStack;
1245 
1246        float maxt, mint;
1247
1248        if (!BoundRay(ray, mint, maxt))
1249                return 0;
1250
1251        Intersectable::NewMail();
1252
1253        Vector3 entp = ray.Extrap(mint);
1254        Vector3 extp = ray.Extrap(maxt);
1255 
1256        VspBspNode *node = mRoot;
1257        VspBspNode *farChild = NULL;
1258       
1259        while (1)
1260        {
1261                if (!node->IsLeaf())
1262                {
1263                        VspBspInterior *in = (VspBspInterior *) node;
1264                       
1265                        Plane3 *splitPlane = in->GetPlane();
1266
1267                        int entSide = splitPlane->Side(entp);
1268                        int extSide = splitPlane->Side(extp);
1269
1270                        Vector3 intersection;
1271
1272                        if (entSide < 0)
1273                        {
1274                                node = in->GetBack();
1275
1276                                if(extSide <= 0) // plane does not split ray => no far child
1277                                        continue;
1278                                       
1279                                farChild = in->GetFront(); // plane splits ray
1280
1281                        } else if (entSide > 0)
1282                        {
1283                                node = in->GetFront();
1284
1285                                if (extSide >= 0) // plane does not split ray => no far child
1286                                        continue;
1287
1288                                farChild = in->GetBack(); // plane splits ray                   
1289                        }
1290                        else // ray and plane are coincident
1291                        {
1292                                // WHAT TO DO IN THIS CASE ?
1293                                //break;
1294                                node = in->GetFront();
1295                                continue;
1296                        }
1297
1298                        // push data for far child
1299                        tStack.push(VspBspRayTraversalData(farChild, extp, maxt));
1300
1301                        // find intersection of ray segment with plane
1302                        float t;
1303                        extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &t);
1304                        maxt *= t;
1305                       
1306                } else // reached leaf => intersection with view cell
1307                {
1308                        VspBspLeaf *leaf = dynamic_cast<VspBspLeaf *>(node);
1309     
1310                        if (!leaf->mViewCell->Mailed())
1311                        {
1312                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
1313                                leaf->mViewCell->Mail();
1314                                ++ hits;
1315                        }
1316                       
1317                        //-- fetch the next far child from the stack
1318                        if (tStack.empty())
1319                                break;
1320     
1321                        entp = extp;
1322                        mint = maxt; // NOTE: need this?
1323
1324                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
1325                                break;
1326
1327                        VspBspRayTraversalData &s = tStack.top();
1328
1329                        node = s.mNode;
1330                        extp = s.mExitPoint;
1331                        maxt = s.mMaxT;
1332
1333                        tStack.pop();
1334                }
1335        }
1336
1337        return hits;
1338}
1339
1340bool VspBspTree::Export(const string filename)
1341{
1342        Exporter *exporter = Exporter::GetExporter(filename);
1343
1344        if (exporter)
1345        {
1346                //exporter->ExportVspBspTree(*this);
1347                return true;
1348        }       
1349
1350        return false;
1351}
1352
1353void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
1354{
1355        stack<VspBspNode *> nodeStack;
1356        nodeStack.push(mRoot);
1357
1358        ViewCell::NewMail();
1359
1360        while (!nodeStack.empty())
1361        {
1362                VspBspNode *node = nodeStack.top();
1363                nodeStack.pop();
1364
1365                if (node->IsLeaf())
1366                {
1367                        ViewCell *viewCell = dynamic_cast<VspBspLeaf *>(node)->mViewCell;
1368
1369                        if (!viewCell->Mailed())
1370                        {
1371                                viewCell->Mail();
1372                                viewCells.push_back(viewCell);
1373                        }
1374                }
1375                else
1376                {
1377                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
1378
1379                        nodeStack.push(interior->mFront);
1380                        nodeStack.push(interior->mBack);
1381                }
1382        }
1383}
1384
1385void VspBspTree::EvaluateViewCellsStats(VspBspViewCellsStatistics &stat) const
1386{
1387        stat.Reset();
1388
1389        stack<VspBspNode *> nodeStack;
1390        nodeStack.push(mRoot);
1391
1392        ViewCell::NewMail();
1393
1394        // exclude root cell
1395        mRootCell->Mail();
1396
1397        while (!nodeStack.empty())
1398        {
1399                VspBspNode *node = nodeStack.top();
1400                nodeStack.pop();
1401
1402                if (node->IsLeaf())
1403                {
1404                        ++ stat.bspLeaves;
1405
1406                        BspViewCell *viewCell = dynamic_cast<VspBspLeaf *>(node)->mViewCell;
1407
1408                        if (!viewCell->Mailed())
1409                        {
1410                                viewCell->Mail();
1411                               
1412                                ++ stat.viewCells;
1413                                const int pvsSize = viewCell->GetPvs().GetSize();
1414
1415                stat.pvs += pvsSize;
1416
1417                                if (pvsSize < 1)
1418                                        ++ stat.emptyPvs;
1419
1420                                if (pvsSize > stat.maxPvs)
1421                                        stat.maxPvs = pvsSize;
1422
1423                                if (pvsSize < stat.minPvs)
1424                                        stat.minPvs = pvsSize;
1425
1426                                //TODO
1427                                //if ((int)viewCell->mVspBspLeaves.size() > stat.maxVspBspLeaves)
1428                                //      stat.maxVspBspLeaves = (int)viewCell->mVspBspLeaves.size();             
1429                        }
1430                }
1431                else
1432                {
1433                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
1434
1435                        nodeStack.push(interior->mFront);
1436                        nodeStack.push(interior->mBack);
1437                }
1438        }
1439}
1440
1441VspBspTreeStatistics &VspBspTree::GetStat()
1442{
1443        return mStat;
1444}
1445
1446float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
1447{
1448        float len = 0;
1449
1450        RayInfoContainer::const_iterator it, it_end = rays.end();
1451
1452        for (it = rays.begin(); it != it_end; ++ it)
1453                len += (*it).SegmentLength();
1454
1455        return len;
1456}
1457
1458int VspBspTree::SplitRays(const Plane3 &plane,
1459                                           RayInfoContainer &rays,
1460                                           RayInfoContainer &frontRays,
1461                                           RayInfoContainer &backRays)
1462{
1463        int splits = 0;
1464
1465        while (!rays.empty())
1466        {
1467                RayInfo bRay = rays.back();
1468                VssRay *ray = bRay.mRay;
1469
1470                rays.pop_back();
1471       
1472                // get classification and receive new t
1473                const int cf =  bRay.ComputeRayIntersection(plane, ray->mT);
1474
1475                switch(cf)
1476                {
1477                case -1:
1478                        backRays.push_back(bRay);
1479                        break;
1480                case 1:
1481                        break;
1482                        frontRays.push_back(bRay);
1483                case 0:
1484                        //-- split ray
1485
1486                        // if start point behind plane
1487                        if (plane.Side(bRay.ExtrapOrigin()) <= 0)
1488                        {
1489                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
1490                                frontRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
1491                        }
1492                        else
1493                        {
1494                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), ray->mT));
1495                                backRays.push_back(RayInfo(ray, ray->mT, bRay.GetMaxT()));
1496                        }
1497                        break;
1498                default:
1499                        Debug << "Should not come here 4" << endl;
1500                        break;
1501                }
1502        }
1503
1504        return splits;
1505}
1506
1507void VspBspTree::ExtractHalfSpaces(VspBspNode *n, vector<Plane3> &halfSpaces) const
1508{
1509        VspBspNode *lastNode;
1510        do
1511        {
1512                lastNode = n;
1513
1514                // want to get planes defining geometry of this node => don't take
1515                // split plane of node itself
1516                n = n->GetParent();
1517               
1518                if (n)
1519                {
1520                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(n);
1521                        Plane3 halfSpace = *dynamic_cast<VspBspInterior *>(interior)->GetPlane();
1522
1523            if (interior->mFront != lastNode)
1524                                halfSpace.ReverseOrientation();
1525
1526                        halfSpaces.push_back(halfSpace);
1527                }
1528        }
1529        while (n);
1530}
1531
1532void VspBspTree::ConstructGeometry(VspBspNode *n, VspBspNodeGeometry &cell) const
1533{
1534        PolygonContainer polys;
1535        ConstructGeometry(n, polys);
1536        cell.mPolys = polys;
1537}
1538
1539void VspBspTree::ConstructGeometry(VspBspNode *n, PolygonContainer &cell) const
1540{
1541        vector<Plane3> halfSpaces;
1542        ExtractHalfSpaces(n, halfSpaces);
1543
1544        PolygonContainer candidatePolys;
1545
1546        // bounded planes are added to the polygons (reverse polygons
1547        // as they have to be outfacing
1548        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
1549        {
1550                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
1551               
1552                if (p->Valid())
1553                {
1554                        candidatePolys.push_back(p->CreateReversePolygon());
1555                        DEL_PTR(p);
1556                }
1557        }
1558
1559        // add faces of bounding box (also could be faces of the cell)
1560        for (int i = 0; i < 6; ++ i)
1561        {
1562                VertexContainer vertices;
1563       
1564                for (int j = 0; j < 4; ++ j)
1565                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
1566
1567                candidatePolys.push_back(new Polygon3(vertices));
1568        }
1569
1570        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
1571        {
1572                // polygon is split by all other planes
1573                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
1574                {
1575                        if (i == j) // polygon and plane are coincident
1576                                continue;
1577
1578                        VertexContainer splitPts;
1579                        Polygon3 *frontPoly, *backPoly;
1580
1581                        const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
1582                       
1583                        switch (cf)
1584                        {
1585                                case Polygon3::SPLIT:
1586                                        frontPoly = new Polygon3();
1587                                        backPoly = new Polygon3();
1588
1589                                        candidatePolys[i]->Split(halfSpaces[j],
1590                                                                                         *frontPoly,
1591                                                                                         *backPoly);
1592
1593                                        DEL_PTR(candidatePolys[i]);
1594
1595                                        if (frontPoly->Valid())
1596                                                candidatePolys[i] = frontPoly;
1597                                        else
1598                                                DEL_PTR(frontPoly);
1599
1600                                        DEL_PTR(backPoly);
1601                                        break;
1602                                case Polygon3::BACK_SIDE:
1603                                        DEL_PTR(candidatePolys[i]);
1604                                        break;
1605                                // just take polygon as it is
1606                                case Polygon3::FRONT_SIDE:
1607                                case Polygon3::COINCIDENT:
1608                                default:
1609                                        break;
1610                        }
1611                }
1612               
1613                if (candidatePolys[i])
1614                        cell.push_back(candidatePolys[i]);
1615        }
1616}
1617
1618void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &cell) const
1619{
1620        /*
1621        TODO
1622        vector<VspBspLeaf *> leaves = vc->mBspLeaves;
1623
1624        vector<VspBspLeaf *>::const_iterator it, it_end = leaves.end();
1625
1626        for (it = leaves.begin(); it != it_end; ++ it)
1627                ConstructGeometry(*it, vcGeom);*/
1628}
1629
1630int VspBspTree::FindNeighbors(VspBspNode *n, vector<VspBspLeaf *> &neighbors,
1631                                                   const bool onlyUnmailed) const
1632{
1633        PolygonContainer cell;
1634
1635        ConstructGeometry(n, cell);
1636
1637        stack<VspBspNode *> nodeStack;
1638        nodeStack.push(mRoot);
1639               
1640        // planes needed to verify that we found neighbor leaf.
1641        vector<Plane3> halfSpaces;
1642        ExtractHalfSpaces(n, halfSpaces);
1643
1644        while (!nodeStack.empty())
1645        {
1646                VspBspNode *node = nodeStack.top();
1647                nodeStack.pop();
1648
1649                if (node->IsLeaf())
1650                {
1651            if (node != n && (!onlyUnmailed || !node->Mailed()))
1652                        {
1653                                // test all planes of current node if candidate really
1654                                // is neighbour
1655                                PolygonContainer neighborCandidate;
1656                                ConstructGeometry(node, neighborCandidate);
1657                               
1658                                bool isAdjacent = true;
1659                                for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
1660                                {
1661                                        const int cf =
1662                                                Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);
1663
1664                                        if (cf == Polygon3::BACK_SIDE)
1665                                                isAdjacent = false;
1666                                }
1667
1668                                if (isAdjacent)
1669                                        neighbors.push_back(dynamic_cast<VspBspLeaf *>(node));
1670
1671                                CLEAR_CONTAINER(neighborCandidate);
1672                        }
1673                }
1674                else
1675                {
1676                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
1677       
1678                        const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);
1679
1680                        if (cf == Polygon3::FRONT_SIDE)
1681                                nodeStack.push(interior->mFront);
1682                        else
1683                                if (cf == Polygon3::BACK_SIDE)
1684                                        nodeStack.push(interior->mBack);
1685                                else
1686                                {
1687                                        // random decision
1688                                        nodeStack.push(interior->mBack);
1689                                        nodeStack.push(interior->mFront);
1690                                }
1691                }
1692        }
1693       
1694        CLEAR_CONTAINER(cell);
1695        return (int)neighbors.size();
1696}
1697
1698VspBspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
1699{
1700    stack<VspBspNode *> nodeStack;
1701        nodeStack.push(mRoot);
1702       
1703        int mask = rand();
1704 
1705        while (!nodeStack.empty())
1706        {
1707                VspBspNode *node = nodeStack.top();
1708                nodeStack.pop();
1709         
1710                if (node->IsLeaf())
1711                {
1712                        return dynamic_cast<VspBspLeaf *>(node);
1713                }
1714                else
1715                {
1716                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
1717                       
1718                        VspBspNode *next;
1719       
1720                        PolygonContainer cell;
1721
1722                        // todo: not very efficient: constructs full cell everytime
1723                        ConstructGeometry(interior, cell);
1724
1725                        const int cf = Polygon3::ClassifyPlane(cell, halfspace);
1726
1727                        if (cf == Polygon3::BACK_SIDE)
1728                                next = interior->mFront;
1729                        else
1730                                if (cf == Polygon3::FRONT_SIDE)
1731                                        next = interior->mFront;
1732                        else
1733                        {
1734                                // random decision
1735                                if (mask & 1)
1736                                        next = interior->mBack;
1737                                else
1738                                        next = interior->mFront;
1739                                mask = mask >> 1;
1740                        }
1741
1742                        nodeStack.push(next);
1743                }
1744        }
1745       
1746        return NULL;
1747}
1748
1749VspBspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
1750{
1751        stack<VspBspNode *> nodeStack;
1752       
1753        nodeStack.push(mRoot);
1754
1755        int mask = rand();
1756       
1757        while (!nodeStack.empty())
1758        {
1759                VspBspNode *node = nodeStack.top();
1760                nodeStack.pop();
1761               
1762                if (node->IsLeaf())
1763                {
1764                        if ( (!onlyUnmailed || !node->Mailed()) )
1765                                return dynamic_cast<VspBspLeaf *>(node);
1766                }
1767                else
1768                {
1769                        VspBspInterior *interior = dynamic_cast<VspBspInterior *>(node);
1770
1771                        // random decision
1772                        if (mask & 1)
1773                                nodeStack.push(interior->mBack);
1774                        else
1775                                nodeStack.push(interior->mFront);
1776
1777                        mask = mask >> 1;
1778                }
1779        }
1780       
1781        return NULL;
1782}
1783
1784int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
1785{
1786        int pvsSize = 0;
1787
1788        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1789
1790        Intersectable::NewMail();
1791
1792        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1793        {
1794                VssRay *ray = (*rit).mRay;
1795               
1796                if (ray->mOriginObject)
1797                {
1798                        if (!ray->mOriginObject->Mailed())
1799                        {
1800                                ray->mOriginObject->Mail();
1801                                ++ pvsSize;
1802                        }
1803                }
1804                if (ray->mTerminationObject)
1805                {
1806                        if (!ray->mTerminationObject->Mailed())
1807                        {
1808                                ray->mTerminationObject->Mail();
1809                                ++ pvsSize;
1810                        }
1811                }
1812        }
1813
1814        return pvsSize;
1815}
1816
1817/*************************************************************
1818 *            VspBspNodeGeometry Implementation                 *
1819 *************************************************************/
1820
1821VspBspNodeGeometry::~VspBspNodeGeometry()
1822{
1823        CLEAR_CONTAINER(mPolys);
1824}
1825
1826float VspBspNodeGeometry::GetArea() const
1827{
1828        return Polygon3::GetArea(mPolys);
1829}
1830
1831void VspBspNodeGeometry::SplitGeometry(VspBspNodeGeometry &front,
1832                                                                        VspBspNodeGeometry &back,
1833                                                                        const VspBspTree &tree,                                         
1834                                                                        const Plane3 &splitPlane) const
1835{       
1836        // get cross section of new polygon
1837        Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);
1838
1839        planePoly = SplitPolygon(planePoly, tree);
1840
1841        //-- plane poly splits all other cell polygons
1842        for (int i = 0; i < (int)mPolys.size(); ++ i)
1843        {
1844                const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
1845                       
1846                // split new polygon with all previous planes
1847                switch (cf)
1848                {
1849                        case Polygon3::SPLIT:
1850                                {
1851                                        Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);
1852
1853                                        Polygon3 *frontPoly = new Polygon3();
1854                                        Polygon3 *backPoly = new Polygon3();
1855                               
1856                                        poly->Split(splitPlane, *frontPoly, *backPoly);
1857
1858                                        DEL_PTR(poly);
1859
1860                                        if (frontPoly->Valid())
1861                                                front.mPolys.push_back(frontPoly);
1862                                        if (backPoly->Valid())
1863                                                back.mPolys.push_back(backPoly);
1864                                }
1865                               
1866                                break;
1867                        case Polygon3::BACK_SIDE:
1868                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));                     
1869                                break;
1870                        case Polygon3::FRONT_SIDE:
1871                                front.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));     
1872                                break;
1873                        case Polygon3::COINCIDENT:
1874                                //front.mPolys.push_back(CreateReversePolygon(mPolys[i]));
1875                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));
1876                                break;
1877                        default:
1878                                break;
1879                }
1880        }
1881
1882        //-- finally add the new polygon to the child cells
1883        if (planePoly)
1884        {
1885                // add polygon with normal pointing into positive half space to back cell
1886                back.mPolys.push_back(planePoly);
1887                // add polygon with reverse orientation to front cell
1888                front.mPolys.push_back(planePoly->CreateReversePolygon());
1889        }
1890
1891        //Debug << "returning new geometry " << mPolys.size() << " f: " << front.mPolys.size() << " b: " << back.mPolys.size() << endl;
1892        //Debug << "old area " << GetArea() << " f: " << front.GetArea() << " b: " << back.GetArea() << endl;
1893}
1894
1895Polygon3 *VspBspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
1896                                                                                const VspBspTree &tree) const
1897{
1898        // polygon is split by all other planes
1899        for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
1900        {
1901                Plane3 plane = mPolys[i]->GetSupportingPlane();
1902
1903                const int cf =
1904                        planePoly->ClassifyPlane(plane, 0.00001f);
1905                       
1906                // split new polygon with all previous planes
1907                switch (cf)
1908                {
1909                        case Polygon3::SPLIT:
1910                                {
1911                                        Polygon3 *frontPoly = new Polygon3();
1912                                        Polygon3 *backPoly = new Polygon3();
1913
1914                                        planePoly->Split(plane, *frontPoly, *backPoly);
1915                                        DEL_PTR(planePoly);
1916
1917                                        if (backPoly->Valid())
1918                                                planePoly = backPoly;
1919                                        else
1920                                                DEL_PTR(backPoly);
1921                                }
1922                                break;
1923                        case Polygon3::FRONT_SIDE:
1924                                DEL_PTR(planePoly);
1925                break;
1926                        // polygon is taken as it is
1927                        case Polygon3::BACK_SIDE:
1928                        case Polygon3::COINCIDENT:
1929                        default:
1930                                break;
1931                }
1932        }
1933
1934        return planePoly;
1935}
1936
1937void VspBspViewCellsStatistics::Print(ostream &app) const
1938{
1939        app << "===== VspBspViewCells statistics ===============\n";
1940
1941        app << setprecision(4);
1942
1943        //app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
1944
1945        app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
1946
1947        app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
1948
1949        app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
1950
1951        app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
1952
1953        app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
1954
1955        app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
1956
1957        app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgVspBspLeaves() << endl;
1958
1959        app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxVspBspLeaves << endl;
1960       
1961        app << "===== END OF VspBspViewCells statistics ==========\n";
1962}
1963
1964BspViewCell *VspBspTree::GetRootCell() const
1965{
1966        return mRootCell;
1967}
Note: See TracBrowser for help on using the repository browser.