source: trunk/VUT/GtpVisibilityPreprocessor/src/ViewCellBsp.cpp @ 423

Revision 423, 61.7 KB checked in by mattausch, 19 years ago (diff)

proceeding with the kd view space partition

Line 
1#include "Plane3.h"
2#include "ViewCellBsp.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
18int BspNode::sMailId = 1;
19
20/** Evaluates split plane classification with respect to the plane's
21        contribution for a minimum number splits in the tree.
22*/
23const float BspTree::sLeastPolySplitsTable[] = {0, 0, 1, 0};
24/** Evaluates split plane classification with respect to the plane's
25        contribution for a balanced tree.
26*/
27const float BspTree::sBalancedPolysTable[] = {1, -1, 0, 0};
28
29/** Evaluates split plane classification with respect to the plane's
30        contribution for a minimum number of ray splits.
31*/
32const float BspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
33/** Evaluates split plane classification with respect to the plane's
34        contribution for balanced rays.
35*/
36const float BspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
37
38int BspTree::sConstructionMethod = FROM_INPUT_VIEW_CELLS;
39
40
41/****************************************************************/
42/*                  class BspNode implementation                */
43/****************************************************************/
44
45BspNode::BspNode():
46mParent(NULL)
47{}
48
49BspNode::BspNode(BspInterior *parent):
50mParent(parent)
51{}
52
53
54bool BspNode::IsRoot() const
55{
56        return mParent == NULL;
57}
58
59BspInterior *BspNode::GetParent()
60{
61        return mParent;
62}
63
64void BspNode::SetParent(BspInterior *parent)
65{
66        mParent = parent;
67}
68
69/****************************************************************/
70/*              class BspInterior implementation                */
71/****************************************************************/
72
73
74BspInterior::BspInterior(const Plane3 &plane):
75mPlane(plane), mFront(NULL), mBack(NULL)
76{}
77
78BspInterior::~BspInterior()
79{
80        DEL_PTR(mFront);
81        DEL_PTR(mBack);
82}
83
84bool BspInterior::IsLeaf() const
85{
86        return false;
87}
88
89BspNode *BspInterior::GetBack()
90{
91        return mBack;
92}
93
94BspNode *BspInterior::GetFront()
95{
96        return mFront;
97}
98
99Plane3 *BspInterior::GetPlane()
100{
101        return &mPlane;
102}
103
104void BspInterior::ReplaceChildLink(BspNode *oldChild, BspNode *newChild)
105{
106        if (mBack == oldChild)
107                mBack = newChild;
108        else
109                mFront = newChild;
110}
111
112void BspInterior::SetupChildLinks(BspNode *b, BspNode *f)
113{
114    mBack = b;
115    mFront = f;
116}
117
118int BspInterior::SplitPolygons(PolygonContainer &polys,
119                                                           PolygonContainer &frontPolys,
120                                                           PolygonContainer &backPolys,
121                                                           PolygonContainer &coincident)
122{
123        Polygon3 *splitPoly = NULL;
124
125        int splits = 0;
126
127#ifdef _Debug
128        Debug << "splitting polygons of node " << this << " with plane " << mPlane << endl;
129#endif
130        while (!polys.empty())
131        {
132                Polygon3 *poly = polys.back();
133                polys.pop_back();
134
135                //Debug << "New polygon with plane: " << poly->GetSupportingPlane() << "\n";
136
137                // classify polygon
138                const int cf = poly->ClassifyPlane(mPlane);
139
140                Polygon3 *front_piece = NULL;
141                Polygon3 *back_piece = NULL;
142       
143                VertexContainer splitVertices;
144
145                switch (cf)
146                {
147                        case Polygon3::COINCIDENT:
148                                coincident.push_back(poly);
149                                break;                 
150                        case Polygon3::FRONT_SIDE:     
151                                frontPolys.push_back(poly);
152                                break;
153                        case Polygon3::BACK_SIDE:
154                                backPolys.push_back(poly);
155                                break;
156                        case Polygon3::SPLIT:
157                                front_piece = new Polygon3(poly->mParent);
158                                back_piece = new Polygon3(poly->mParent);
159
160                                //-- split polygon into front and back part
161                                poly->Split(mPlane,
162                                                        *front_piece,
163                                                        *back_piece,
164                                                        splitVertices);
165                                       
166                                ++ splits; // increase number of splits
167
168                                //-- inherit rays from parent polygon for blocked ray criterium
169                                poly->InheritRays(*front_piece, *back_piece);
170                                //Debug << "p: " << poly->mPiercingRays.size() << " f: " << front_piece->mPiercingRays.size() << " b: " << back_piece->mPiercingRays.size() << endl;
171                               
172                                // check if polygons still valid
173                                if (front_piece->Valid())
174                                        frontPolys.push_back(front_piece);
175                                else
176                                        DEL_PTR(front_piece);
177                               
178                                if (back_piece->Valid())
179                                        backPolys.push_back(back_piece);
180                                else                           
181                                        DEL_PTR(back_piece);
182                               
183#ifdef _DEBUG
184                                Debug << "split " << *poly << endl << *front_piece << endl << *back_piece << endl;
185#endif
186                                DEL_PTR(poly);                 
187                                break;
188                        default:
189                Debug << "SHOULD NEVER COME HERE\n";
190                                break;
191                }
192        }
193
194        return splits;
195}
196
197/****************************************************************/
198/*                  class BspLeaf implementation                */
199/****************************************************************/
200BspLeaf::BspLeaf(): mViewCell(NULL)
201{
202}
203
204BspLeaf::BspLeaf(BspViewCell *viewCell):
205mViewCell(viewCell)
206{
207}
208
209BspLeaf::BspLeaf(BspInterior *parent):
210BspNode(parent), mViewCell(NULL)
211{}
212
213BspLeaf::BspLeaf(BspInterior *parent, BspViewCell *viewCell):
214BspNode(parent), mViewCell(viewCell)
215{
216}
217
218BspViewCell *BspLeaf::GetViewCell() const
219{
220        return mViewCell;
221}
222
223void BspLeaf::SetViewCell(BspViewCell *viewCell)
224{
225        mViewCell = viewCell;
226}
227
228bool BspLeaf::IsLeaf() const
229{
230        return true;
231}
232
233void BspLeaf::AddToPvs(const BoundedRayContainer &rays,
234                                           int &sampleContributions,
235                                           int &contributingSamples,
236                                           bool storeLeavesWithRays)
237{
238        sampleContributions = 0;
239        contributingSamples = 0;
240
241    BoundedRayContainer::const_iterator it, it_end = rays.end();
242
243        // add contributions from samples to the PVS
244        for (it = rays.begin(); it != it_end; ++ it)
245        {
246                int contribution = 0;
247                Ray *ray = (*it)->mRay;
248                       
249                if (!ray->intersections.empty())
250                        contribution += mViewCell->GetPvs().AddSample(ray->intersections[0].mObject);
251               
252                if (ray->sourceObject.mObject)
253                        contribution += mViewCell->GetPvs().AddSample(ray->sourceObject.mObject);
254
255                if (contribution > 0)
256                {
257                        sampleContributions += contribution;
258                        ++ contributingSamples;
259                }
260               
261                if (storeLeavesWithRays)
262                        // warning: intersections not ordered
263                        ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
264        }
265}
266
267
268/****************************************************************/
269/*                  class BspTree implementation                */
270/****************************************************************/
271
272BspTree::BspTree(BspViewCell *viewCell):
273mRootCell(viewCell),
274mRoot(NULL),
275mGenerateViewCells(false),
276//-- factors for bsp tree split plane heuristics
277mVerticalSplitsFactor(1.0f),
278mLargestPolyAreaFactor(1.0f),
279mBlockedRaysFactor(1.0f),
280mLeastRaySplitsFactor(1.0f),
281mBalancedRaysFactor(1.0f),
282mPvsFactor(1.0f),
283mLeastSplitsFactor(1.0f),
284mBalancedPolysFactor(1.0f),
285mBalancedViewCellsFactor(1.0f),
286//------------------------------------
287mStoreLeavesWithRays(false),
288mPvsUseArea(true)
289{
290        Randomize(); // initialise random generator for heuristics
291
292        //-- termination criteria for autopartition
293        environment->GetIntValue("BspTree.Termination.maxDepth", mTermMaxDepth);
294        environment->GetIntValue("BspTree.Termination.minPvs", mTermMinPvs);
295        environment->GetIntValue("BspTree.Termination.maxPolygons", mTermMaxPolygons);
296        environment->GetIntValue("BspTree.Termination.maxRays", mTermMaxRays);
297        environment->GetFloatValue("BspTree.Termination.minArea", mTermMinArea);       
298        environment->GetFloatValue("BspTree.Termination.maxRayContribution", mTermMaxRayContribution);
299        environment->GetFloatValue("BspTree.Termination.maxAccRayLenght", mTermMaxAccRayLength);
300
301        //-- termination criteria for axis aligned split
302        environment->GetFloatValue("BspTree.Termination.AxisAligned.ct_div_ci", mCtDivCi);
303        environment->GetFloatValue("BspTree.Termination.AxisAligned.maxCostRatio", mMaxCostRatio);
304        environment->GetIntValue("BspTree.Termination.AxisAligned.maxPolys",
305                                                         mTermMaxPolysForAxisAligned);
306        environment->GetIntValue("BspTree.Termination.AxisAligned.maxRays",
307                                                         mTermMaxRaysForAxisAligned);
308        environment->GetIntValue("BspTree.Termination.AxisAligned.maxObjects",
309                                                         mTermMaxObjectsForAxisAligned);
310        //-- partition criteria
311        environment->GetIntValue("BspTree.maxPolyCandidates", mMaxPolyCandidates);
312        environment->GetIntValue("BspTree.maxRayCandidates", mMaxRayCandidates);
313        environment->GetIntValue("BspTree.splitPlaneStrategy", mSplitPlaneStrategy);
314        environment->GetFloatValue("BspTree.AxisAligned.splitBorder", mSplitBorder);
315
316        environment->GetFloatValue("BspTree.Construction.sideTolerance", Vector3::sDistTolerance);
317        Vector3::sDistToleranceSqrt = Vector3::sDistTolerance * Vector3::sDistTolerance;
318
319        // post processing stuff
320        environment->GetIntValue("ViewCells.PostProcessing.minPvsDif", mMinPvsDif);
321        environment->GetIntValue("ViewCells.PostProcessing.minPvs", mMinPvs);
322        environment->GetIntValue("ViewCells.PostProcessing.maxPvs", mMaxPvs);
323
324    Debug << "BSP max depth: " << mTermMaxDepth << endl;
325        Debug << "BSP min PVS: " << mTermMinPvs << endl;
326        Debug << "BSP min area: " << mTermMinArea << endl;
327        Debug << "BSP max polys: " << mTermMaxPolygons << endl;
328        Debug << "BSP max rays: " << mTermMaxRays << endl;
329        Debug << "BSP max polygon candidates: " << mMaxPolyCandidates << endl;
330        Debug << "BSP max plane candidates: " << mMaxRayCandidates << endl;
331
332        Debug << "Split plane strategy: ";
333        if (mSplitPlaneStrategy & RANDOM_POLYGON)
334                Debug << "random polygon ";
335        if (mSplitPlaneStrategy & AXIS_ALIGNED)
336                Debug << "axis aligned ";
337        if (mSplitPlaneStrategy & LEAST_SPLITS)     
338                Debug << "least splits ";
339        if (mSplitPlaneStrategy & BALANCED_POLYS)
340                Debug << "balanced polygons ";
341        if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
342                Debug << "balanced view cells ";
343        if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
344                Debug << "largest polygon area ";
345        if (mSplitPlaneStrategy & VERTICAL_AXIS)
346                Debug << "vertical axis ";
347        if (mSplitPlaneStrategy & BLOCKED_RAYS)
348                Debug << "blocked rays ";
349        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
350                Debug << "least ray splits ";
351        if (mSplitPlaneStrategy & BALANCED_RAYS)
352                Debug << "balanced rays ";
353        if (mSplitPlaneStrategy & PVS)
354                Debug << "pvs";
355
356        Debug << endl;
357}
358
359void BspTree::ParseEnvironment()
360{
361        //-- parse bsp cell tree construction method
362        char constructionMethodStr[60];
363       
364        environment->GetStringValue("BspTree.Construction.input", constructionMethodStr);
365
366        sConstructionMethod = FROM_INPUT_VIEW_CELLS;
367       
368        if (strcmp(constructionMethodStr, "fromViewCells") == 0)
369                sConstructionMethod = FROM_INPUT_VIEW_CELLS;
370        else if (strcmp(constructionMethodStr, "fromSceneGeometry") == 0)
371                sConstructionMethod = FROM_SCENE_GEOMETRY;
372        else if (strcmp(constructionMethodStr, "fromSamples") == 0)
373                sConstructionMethod = FROM_SAMPLES;
374        else
375        {
376                cerr << "Wrong construction method " << constructionMethodStr << endl;
377                exit(1);
378    }
379
380        Debug << "Construction method: " << constructionMethodStr << endl;
381}
382
383const BspTreeStatistics &BspTree::GetStatistics() const
384{
385        return mStat;
386}
387
388void BspTreeStatistics::Print(ostream &app) const
389{
390        app << "===== BspTree statistics ===============\n";
391
392        app << setprecision(4);
393
394        app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
395
396        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
397
398        app << "#N_INTERIORS ( Number of interior nodes )\n" << Interior() << "\n";
399
400        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
401
402        app << "#N_SPLITS ( Number of splits )\n" << splits << "\n";
403
404        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
405        maxDepthNodes * 100 / (double)Leaves() << endl;
406
407        app << "#N_PMAXDEPTH ( Maximal reached depth )\n" << maxDepth << endl;
408
409        app << "#N_PMINDEPTH ( Minimal reached depth )\n" << minDepth << endl;
410
411        app << "#AVGDEPTH ( average depth )\n" << AvgDepth() << endl;
412
413        app << "#N_INPUT_POLYGONS (number of input polygons )\n" <<     polys << endl;
414
415        //app << "#N_PVS: " << pvs << endl;
416
417        app << "#N_ROUTPUT_INPUT_POLYGONS ( ratio polygons after subdivision / input polygons )\n" <<
418                 (polys + splits) / (double)polys << endl;
419       
420        app << "===== END OF BspTree statistics ==========\n";
421}
422
423
424BspTree::~BspTree()
425{
426        DEL_PTR(mRoot);
427}
428
429void BspTree::InsertViewCell(ViewCell *viewCell)
430{
431        PolygonContainer *polys = new PolygonContainer();
432
433        // extract polygons that guide the split process
434        mStat.polys += AddMeshToPolygons(viewCell->GetMesh(), *polys, viewCell);
435        mBox.Include(viewCell->GetBox()); // add to BSP aabb
436
437        InsertPolygons(polys);
438}
439
440void BspTree::InsertPolygons(PolygonContainer *polys)
441{       
442        std::stack<BspTraversalData> tStack;
443       
444        // traverse existing tree or create new tree
445    if (!mRoot)
446                mRoot = new BspLeaf();
447
448        tStack.push(BspTraversalData(mRoot, polys, 0, mRootCell, new BoundedRayContainer(), 0,
449                                                                 mBox.SurfaceArea(), new BspNodeGeometry()));
450
451        while (!tStack.empty())
452        {
453                // filter polygons donw the tree
454                BspTraversalData tData = tStack.top();
455            tStack.pop();
456                       
457                if (!tData.mNode->IsLeaf())
458                {
459                        BspInterior *interior = dynamic_cast<BspInterior *>(tData.mNode);
460
461                        //-- filter view cell polygons down the tree until a leaf is reached
462                        if (!tData.mPolygons->empty())
463                        {
464                                PolygonContainer *frontPolys = new PolygonContainer();
465                                PolygonContainer *backPolys = new PolygonContainer();
466                                PolygonContainer coincident;
467
468                                int splits = 0;
469               
470                                // split viewcell polygons with respect to split plane
471                                splits += interior->SplitPolygons(*tData.mPolygons,
472                                                                                                  *frontPolys,
473                                                                                                  *backPolys,
474                                                                                                  coincident);
475                               
476                                // extract view cells associated with the split polygons
477                                ViewCell *frontViewCell = mRootCell;
478                                ViewCell *backViewCell = mRootCell;
479                       
480                                BspTraversalData frontData(interior->GetFront(),
481                                                                                   frontPolys,
482                                                                                   tData.mDepth + 1,
483                                                                                   mRootCell,   
484                                                                                   tData.mRays,
485                                                                                   tData.mPvs,
486                                                                                   mBox.SurfaceArea(),
487                                                                                   new BspNodeGeometry());
488
489                                BspTraversalData backData(interior->GetBack(),
490                                                                                  backPolys,
491                                                                                  tData.mDepth + 1,
492                                                                                  mRootCell,   
493                                                                                  tData.mRays,
494                                                                                  tData.mPvs,
495                                                                                  mBox.SurfaceArea(),
496                                                                                  new BspNodeGeometry());
497
498                                if (!mGenerateViewCells)
499                                {
500                                        ExtractViewCells(frontData,
501                                                                         backData,
502                                                                         coincident,
503                                                                         interior->mPlane);
504                                }
505
506                                // don't need coincident polygons anymore
507                                CLEAR_CONTAINER(coincident);
508
509                                mStat.splits += splits;
510
511                                // push the children on the stack
512                                tStack.push(frontData);
513                                tStack.push(backData);
514                        }
515
516                        // cleanup
517                        DEL_PTR(tData.mPolygons);
518                }
519                else
520                {
521                        // reached leaf => subdivide current viewcell
522                        BspNode *subRoot = Subdivide(tStack, tData);
523                }
524        }
525}
526
527int BspTree::AddMeshToPolygons(Mesh *mesh, PolygonContainer &polys, MeshInstance *parent)
528{
529        FaceContainer::const_iterator fi;
530       
531        // copy the face data to polygons
532        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
533        {
534                Polygon3 *poly = new Polygon3((*fi), mesh);
535               
536                if (poly->Valid())
537                {
538                        poly->mParent = parent; // set parent intersectable
539                        polys.push_back(poly);
540                }
541                else
542                        DEL_PTR(poly);
543        }
544        return (int)mesh->mFaces.size();
545}
546
547int BspTree::AddToPolygonSoup(const ViewCellContainer &viewCells,
548                                                          PolygonContainer &polys,
549                                                          int maxObjects)
550{
551        int limit = (maxObjects > 0) ?
552                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
553 
554        int polysSize = 0;
555
556        for (int i = 0; i < limit; ++ i)
557        {
558                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
559                {
560                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
561                        polysSize += AddMeshToPolygons(viewCells[i]->GetMesh(), polys, viewCells[i]);
562                }
563        }
564
565        return polysSize;
566}
567
568int BspTree::AddToPolygonSoup(const ObjectContainer &objects, PolygonContainer &polys, int maxObjects)
569{
570        int limit = (maxObjects > 0) ? Min((int)objects.size(), maxObjects) : (int)objects.size();
571 
572        for (int i = 0; i < limit; ++i)
573        {
574                Intersectable *object = objects[i];//*it;
575                Mesh *mesh = NULL;
576
577                switch (object->Type()) // extract the meshes
578                {
579                case Intersectable::MESH_INSTANCE:
580                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
581                        break;
582                case Intersectable::VIEW_CELL:
583                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
584                        break;
585                        // TODO: handle transformed mesh instances
586                default:
587                        Debug << "intersectable type not supported" << endl;
588                        break;
589                }
590               
591        if (mesh) // copy the mesh data to polygons
592                {
593                        mBox.Include(object->GetBox()); // add to BSP tree aabb
594                        AddMeshToPolygons(mesh, polys, mRootCell);
595                }
596        }
597
598        return (int)polys.size();
599}
600
601void BspTree::Construct(const ViewCellContainer &viewCells)
602{
603        mStat.nodes = 1;
604        mBox.Initialize();      // initialise bsp tree bounding box
605
606        // copy view cell meshes into one big polygon soup
607        PolygonContainer *polys = new PolygonContainer();
608        mStat.polys = AddToPolygonSoup(viewCells, *polys);
609
610        // construct tree from the view cell polygons
611        Construct(polys, new BoundedRayContainer());
612}
613
614
615void BspTree::Construct(const ObjectContainer &objects)
616{
617        mStat.nodes = 1;
618        mBox.Initialize();      // initialise bsp tree bounding box
619       
620        PolygonContainer *polys = new PolygonContainer();
621       
622        // copy mesh instance polygons into one big polygon soup
623        mStat.polys = AddToPolygonSoup(objects, *polys);
624
625        // construct tree from polygon soup
626        Construct(polys, new BoundedRayContainer());
627}
628
629void BspTree::Construct(const RayContainer &sampleRays)
630{
631    mStat.nodes = 1;
632        mBox.Initialize();      // initialise BSP tree bounding box
633       
634        PolygonContainer *polys = new PolygonContainer();
635        BoundedRayContainer *rays = new BoundedRayContainer();
636
637        RayContainer::const_iterator rit, rit_end = sampleRays.end();
638
639        long startTime = GetTime();
640
641        Debug << "**** Extracting polygons from rays ****\n";
642
643        std::map<Face *, Polygon3 *> facePolyMap;
644
645        //-- extract polygons intersected by the rays
646        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
647        {
648                Ray *ray = *rit;
649       
650                // get ray-face intersection. Store polygon representing the rays together
651                // with rays intersecting the face.
652                if (!ray->intersections.empty())
653                {
654                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->intersections[0].mObject);
655                        Face *face = obj->GetMesh()->mFaces[ray->intersections[0].mFace];
656
657                        std::map<Face *, Polygon3 *>::iterator it = facePolyMap.find(face);
658
659                        if (it != facePolyMap.end())
660                        {
661                                //store rays if needed for heuristics
662                                if (mSplitPlaneStrategy & BLOCKED_RAYS)
663                                        (*it).second->mPiercingRays.push_back(ray);
664                        }
665                        else
666                        {       //store rays if needed for heuristics
667                                Polygon3 *poly = new Polygon3(face, obj->GetMesh());
668                                poly->mParent = obj;
669                                polys->push_back(poly);
670
671                                if (mSplitPlaneStrategy & BLOCKED_RAYS)
672                                        poly->mPiercingRays.push_back(ray);
673
674                                facePolyMap[face] = poly;
675                        }
676                }
677        }
678       
679        facePolyMap.clear();
680
681        // compue bounding box
682        Polygon3::IncludeInBox(*polys, mBox);
683
684        //-- store rays
685        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
686        {
687                Ray *ray = *rit;
688        ray->SetId(-1); // reset id
689
690                float minT, maxT;
691                if (BoundRay(*ray, minT, maxT))
692                        rays->push_back(new BoundedRay(ray, minT, maxT));
693        }
694
695        mStat.polys = (int)polys->size();
696
697        Debug << "**** Finished polygon extraction ****" << endl;
698        Debug << (int)polys->size() << " polys extracted from " << (int)sampleRays.size() << " rays" << endl;
699        Debug << "extraction time: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
700
701        Construct(polys, rays);
702}
703
704void BspTree::Construct(PolygonContainer *polys, BoundedRayContainer *rays)
705{
706        std::stack<BspTraversalData> tStack;
707
708        mRoot = new BspLeaf();
709
710        BspNodeGeometry *cell = new BspNodeGeometry();
711        ConstructGeometry(mRoot, *cell);
712
713        BspTraversalData tData(mRoot, polys, 0, mRootCell, rays,
714                                                   ComputePvsSize(*rays), cell->GetArea(), cell);
715
716        tStack.push(tData);
717
718        mStat.Start();
719        cout << "**** Contructing bsp tree ****\n";
720
721        while (!tStack.empty())
722        {
723                tData = tStack.top();
724            tStack.pop();
725
726                // subdivide leaf node
727                BspNode *subRoot = Subdivide(tStack, tData);
728        }
729
730        mStat.Stop();
731}
732
733bool BspTree::TerminationCriteriaMet(const BspTraversalData &data) const
734{
735        return
736                (((int)data.mPolygons->size() <= mTermMaxPolygons) ||
737                 ((int)data.mRays->size() <= mTermMaxRays) ||
738                 (data.mPvs <= mTermMinPvs) ||
739                 (data.mArea <= mTermMinArea) ||
740                 (data.mDepth >= mTermMaxDepth) ||
741                 (data.GetAvgRayContribution() < mTermMaxRayContribution));
742}
743
744BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
745{
746        //-- terminate traversal 
747        if (TerminationCriteriaMet(tData))             
748        {
749                BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
750       
751                BspViewCell *viewCell;
752
753                // generate new view cell for each leaf
754                if (mGenerateViewCells)
755                {
756                        viewCell = dynamic_cast<BspViewCell *>(ViewCell::Generate());
757                }
758                else
759                {
760                        // add view cell to leaf
761                        viewCell = dynamic_cast<BspViewCell *>(tData.mViewCell);
762                }
763
764                leaf->SetViewCell(viewCell);
765                viewCell->mBspLeaves.push_back(leaf);
766
767                //-- add pvs
768                if (viewCell != mRootCell)
769                {
770                        int conSamp = 0, sampCon = 0;
771                        leaf->AddToPvs(*tData.mRays, conSamp, sampCon);
772                       
773                        mStat.contributingSamples += conSamp;
774                        mStat.sampleContributions += sampCon;
775                }
776
777                EvaluateLeafStats(tData);
778               
779                //-- clean up
780               
781                // discard polygons
782                CLEAR_CONTAINER(*tData.mPolygons);
783                // discard rays
784                CLEAR_CONTAINER(*tData.mRays);
785
786                DEL_PTR(tData.mPolygons);
787                DEL_PTR(tData.mRays);
788                DEL_PTR(tData.mGeometry);
789               
790                return leaf;
791        }
792
793        //-- continue subdivision
794        PolygonContainer coincident;
795       
796        PolygonContainer *frontPolys = new PolygonContainer();
797        PolygonContainer *backPolys = new PolygonContainer();
798
799        BoundedRayContainer *frontRays = new BoundedRayContainer();
800        BoundedRayContainer *backRays = new BoundedRayContainer();
801       
802        BspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell,
803                                                                new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
804        BspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mRootCell,
805                                                           new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
806
807        // create new interior node and two leaf nodes
808        BspInterior *interior = SubdivideNode(tData,
809                                                                                  tFrontData,
810                                                                                  tBackData,
811                                                                                  coincident);
812
813#ifdef _DEBUG   
814        if (frontPolys->empty() && backPolys->empty() && (coincident.size() > 2))
815        {       for (PolygonContainer::iterator it = coincident.begin(); it != coincident.end(); ++it)
816                        Debug << (*it) << " " << (*it)->GetArea() << " " << (*it)->mParent << endl ;
817                Debug << endl;}
818#endif
819
820        // extract view cells from coincident polygons according to plane normal
821    // only if front or back polygons are empty
822        if (!mGenerateViewCells)
823        {
824                ExtractViewCells(tFrontData,
825                                                 tBackData,
826                                                 coincident,
827                                                 interior->mPlane);                     
828        }
829
830        // don't need coincident polygons anymory
831        CLEAR_CONTAINER(coincident);
832
833        // push the children on the stack
834        tStack.push(tFrontData);
835        tStack.push(tBackData);
836
837        // cleanup
838        DEL_PTR(tData.mNode);
839        DEL_PTR(tData.mPolygons);
840        DEL_PTR(tData.mRays);
841        DEL_PTR(tData.mGeometry);               
842
843        return interior;
844}
845
846void BspTree::ExtractViewCells(BspTraversalData &frontData,
847                                                           BspTraversalData &backData,
848                                                           const PolygonContainer &coincident,
849                                                           const Plane3 splitPlane) const
850{
851        // if not empty, tree is further subdivided => don't have to find view cell
852        bool foundFront = !frontData.mPolygons->empty();
853        bool foundBack = !frontData.mPolygons->empty();
854
855        PolygonContainer::const_iterator it =
856                coincident.begin(), it_end = coincident.end();
857
858        //-- find first view cells in front and back leafs
859        for (; !(foundFront && foundBack) && (it != it_end); ++ it)
860        {
861                if (DotProd((*it)->GetNormal(), splitPlane.mNormal) > 0)
862                {
863                        backData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
864                        foundBack = true;
865                }
866                else
867                {
868                        frontData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
869                        foundFront = true;
870                }
871        }
872}
873
874BspInterior *BspTree::SubdivideNode(BspTraversalData &tData,
875                                                                        BspTraversalData &frontData,
876                                                                        BspTraversalData &backData,
877                                                                        PolygonContainer &coincident)
878{
879        mStat.nodes += 2;
880       
881        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
882        // select subdivision plane
883        BspInterior *interior =
884                new BspInterior(SelectPlane(leaf, tData));
885
886#ifdef _DEBUG
887        Debug << interior << endl;
888#endif
889       
890        // subdivide rays into front and back rays
891        SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);
892       
893        // subdivide polygons with plane
894        mStat.splits += interior->SplitPolygons(*tData.mPolygons,
895                                                    *frontData.mPolygons,
896                                                                                        *backData.mPolygons,
897                                                                                        coincident);
898
899    // compute pvs
900        frontData.mPvs = ComputePvsSize(*frontData.mRays);
901        backData.mPvs = ComputePvsSize(*backData.mRays);
902
903        // split geometry and compute area
904        if (1)
905        {
906                tData.mGeometry->SplitGeometry(*frontData.mGeometry,
907                                                                           *backData.mGeometry,
908                                                                           *this,
909                                                                           interior->mPlane);
910       
911               
912                frontData.mArea = frontData.mGeometry->GetArea();
913                backData.mArea = backData.mGeometry->GetArea();
914        }
915
916        // compute accumulated ray length
917        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
918        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
919
920        //-- create front and back leaf
921
922        BspInterior *parent = leaf->GetParent();
923
924        // replace a link from node's parent
925        if (!leaf->IsRoot())
926        {
927                parent->ReplaceChildLink(leaf, interior);
928                interior->SetParent(parent);
929        }
930        else // new root
931        {
932                mRoot = interior;
933        }
934
935        // and setup child links
936        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
937       
938        frontData.mNode = interior->mFront;
939        backData.mNode = interior->mBack;
940       
941        //DEL_PTR(leaf);
942        return interior;
943}
944
945void BspTree::SortSplitCandidates(const PolygonContainer &polys,
946                                                                  const int axis,
947                                                                  vector<SortableEntry> &splitCandidates) const
948{
949  splitCandidates.clear();
950 
951  int requestedSize = 2 * (int)polys.size();
952  // creates a sorted split candidates array 
953  splitCandidates.reserve(requestedSize);
954 
955  PolygonContainer::const_iterator it, it_end = polys.end();
956
957  AxisAlignedBox3 box;
958
959  // insert all queries
960  for(it = polys.begin(); it != it_end; ++ it)
961  {
962          box.Initialize();
963          box.Include(*(*it));
964         
965          splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
966      splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
967  }
968 
969  stable_sort(splitCandidates.begin(), splitCandidates.end());
970}
971
972
973float BspTree::BestCostRatio(const PolygonContainer &polys,
974                                                         const AxisAlignedBox3 &box,
975                                                         const int axis,
976                                                         float &position,
977                                                         int &objectsBack,
978                                                         int &objectsFront) const
979{
980        vector<SortableEntry> splitCandidates;
981
982        SortSplitCandidates(polys, axis, splitCandidates);
983       
984        // go through the lists, count the number of objects left and right
985        // and evaluate the following cost funcion:
986        // C = ct_div_ci  + (ol + or)/queries
987       
988        int objectsLeft = 0, objectsRight = (int)polys.size();
989       
990        float minBox = box.Min(axis);
991        float maxBox = box.Max(axis);
992        float boxArea = box.SurfaceArea();
993 
994        float minBand = minBox + mSplitBorder * (maxBox - minBox);
995        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
996       
997        float minSum = 1e20f;
998        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
999
1000        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci)
1001        {
1002                switch ((*ci).type)
1003                {
1004                        case SortableEntry::POLY_MIN:
1005                                ++ objectsLeft;
1006                                break;
1007                        case SortableEntry::POLY_MAX:
1008                            -- objectsRight;
1009                                break;
1010                        default:
1011                                break;
1012                }
1013               
1014                if ((*ci).value > minBand && (*ci).value < maxBand)
1015                {
1016                        AxisAlignedBox3 lbox = box;
1017                        AxisAlignedBox3 rbox = box;
1018                        lbox.SetMax(axis, (*ci).value);
1019                        rbox.SetMin(axis, (*ci).value);
1020
1021                        float sum = objectsLeft * lbox.SurfaceArea() +
1022                                                objectsRight * rbox.SurfaceArea();
1023     
1024                        if (sum < minSum)
1025                        {
1026                                minSum = sum;
1027                                position = (*ci).value;
1028
1029                                objectsBack = objectsLeft;
1030                                objectsFront = objectsRight;
1031                        }
1032                }
1033        }
1034 
1035        float oldCost = (float)polys.size();
1036        float newCost = mCtDivCi + minSum / boxArea;
1037        float ratio = newCost / oldCost;
1038
1039
1040#if 0
1041  Debug << "====================" << endl;
1042  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
1043      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
1044#endif
1045  return ratio;
1046}
1047
1048bool BspTree::SelectAxisAlignedPlane(Plane3 &plane,
1049                                                                         const PolygonContainer &polys) const
1050{
1051        AxisAlignedBox3 box;
1052        box.Initialize();
1053       
1054        // create bounding box of region
1055        Polygon3::IncludeInBox(polys, box);
1056       
1057        int objectsBack = 0, objectsFront = 0;
1058        int axis = 0;
1059        float costRatio = MAX_FLOAT;
1060        Vector3 position;
1061
1062        //-- area subdivision
1063        for (int i = 0; i < 3; ++ i)
1064        {
1065                float p = 0;
1066                float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
1067               
1068                if (r < costRatio)
1069                {
1070                        costRatio = r;
1071                        axis = i;
1072                        position = p;
1073                }
1074        }
1075       
1076        if (costRatio >= mMaxCostRatio)
1077                return false;
1078
1079        Vector3 norm(0,0,0); norm[axis] = 1.0f;
1080        plane = Plane3(norm, position);
1081
1082        return true;
1083}
1084
1085Plane3 BspTree::SelectPlane(BspLeaf *leaf, BspTraversalData &data)
1086{
1087        if (data.mPolygons->empty() && data.mRays->empty())
1088        {
1089                Debug << "Warning: No autopartition polygon candidate available\n";
1090       
1091                // return axis aligned split
1092                AxisAlignedBox3 box;
1093                box.Initialize();
1094       
1095                // create bounding box of region
1096                Polygon3::IncludeInBox(*data.mPolygons, box);
1097
1098                const int axis = box.Size().DrivingAxis();
1099                const Vector3 position = (box.Min()[axis] + box.Max()[axis])*0.5f;
1100
1101                Vector3 norm(0,0,0); norm[axis] = 1.0f;
1102                return Plane3(norm, position);
1103        }
1104       
1105        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
1106                ((int)data.mPolygons->size() > mTermMaxPolysForAxisAligned) &&
1107                ((int)data.mRays->size() > mTermMaxRaysForAxisAligned) &&
1108                ((mTermMaxObjectsForAxisAligned < 0) ||
1109                  (Polygon3::ParentObjectsSize(*data.mPolygons) > mTermMaxObjectsForAxisAligned)))
1110        {
1111                Plane3 plane;
1112                if (SelectAxisAlignedPlane(plane, *data.mPolygons))
1113                        return plane;
1114        }
1115
1116        // simplest strategy: just take next polygon
1117        if (mSplitPlaneStrategy & RANDOM_POLYGON)
1118        {
1119        if (!data.mPolygons->empty())
1120                {
1121                        Polygon3 *nextPoly = (*data.mPolygons)[Random((int)data.mPolygons->size())];
1122                        return nextPoly->GetSupportingPlane();
1123                }
1124                else
1125                {
1126                        const int candidateIdx = Random((int)data.mRays->size());
1127                        BoundedRay *bRay = (*data.mRays)[candidateIdx];
1128
1129                        Ray *ray = bRay->mRay;
1130                                               
1131                        const Vector3 minPt = ray->Extrap(bRay->mMinT);
1132                        const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
1133
1134                        const Vector3 pt = (maxPt + minPt) * 0.5;
1135
1136                        const Vector3 normal = ray->GetDir();
1137                       
1138                        return Plane3(normal, pt);
1139                }
1140
1141                return Plane3();
1142        }
1143
1144        // use heuristics to find appropriate plane
1145        return SelectPlaneHeuristics(leaf, data);
1146}
1147
1148Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf, BspTraversalData &data)
1149{
1150        float lowestCost = MAX_FLOAT;
1151        Plane3 bestPlane;
1152        Plane3 plane;
1153
1154        int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
1155       
1156        int candidateIdx = limit;
1157
1158        for (int i = 0; i < limit; ++ i)
1159        {
1160                candidateIdx = GetNextCandidateIdx(candidateIdx, *data.mPolygons);
1161               
1162                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
1163                // evaluate current candidate
1164                const float candidateCost =
1165                        SplitPlaneCost(poly->GetSupportingPlane(), data);
1166
1167                if (candidateCost < lowestCost)
1168                {
1169                        bestPlane = poly->GetSupportingPlane();
1170                        lowestCost = candidateCost;
1171                }
1172        }
1173       
1174        //Debug << "lowest: " << lowestCost << endl;
1175
1176        //-- choose candidate planes extracted from rays
1177        // we currently use two methods
1178        // 1) take 3 ray endpoints, where two are minimum and one a maximum
1179        //    point or the other way round
1180        // 2) take plane normal as plane normal and the midpoint of the ray.
1181        //    PROBLEM: does not resemble any point where visibility is likely to change
1182        const BoundedRayContainer *rays = data.mRays;
1183
1184        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
1185        {
1186                candidateIdx = Random((int)rays->size());
1187                BoundedRay *bRay = (*rays)[candidateIdx];
1188
1189                Ray *ray = bRay->mRay;
1190                                               
1191                const Vector3 minPt = ray->Extrap(bRay->mMinT);
1192                const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
1193
1194                const Vector3 pt = (maxPt + minPt) * 0.5;
1195
1196                const Vector3 normal = ray->GetDir();
1197                       
1198                plane = Plane3(normal, pt);
1199       
1200                const float candidateCost = SplitPlaneCost(plane, data);
1201
1202                if (candidateCost < lowestCost)
1203                {
1204                        bestPlane = plane;
1205                       
1206                        lowestCost = candidateCost;
1207                }
1208        }
1209
1210        //Debug << "lowest: " << lowestCost << endl;
1211        for (int i = 0; i < mMaxRayCandidates / 2; ++ i)
1212        {
1213                Vector3 pt[3];
1214                int idx[3];
1215                int cmaxT = 0;
1216                int cminT = 0;
1217                bool chooseMin = false;
1218
1219                for (int j = 0; j < 3; j ++)
1220                {
1221                        idx[j] = Random((int)rays->size() * 2);
1222                               
1223                        if (idx[j] >= (int)rays->size())
1224                        {
1225                                idx[j] -= (int)rays->size();
1226                               
1227                                chooseMin = (cminT < 2);
1228                        }
1229                        else
1230                                chooseMin = (cmaxT < 2);
1231
1232                        BoundedRay *bRay = (*rays)[idx[j]];
1233                        pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) : bRay->mRay->Extrap(bRay->mMaxT);
1234                }       
1235                       
1236                plane = Plane3(pt[0], pt[1], pt[2]);
1237
1238                const float candidateCost = SplitPlaneCost(plane, data);
1239
1240                if (candidateCost < lowestCost)
1241                {
1242                        //Debug << "choose ray plane 2: " << candidateCost << endl;
1243                        bestPlane = plane;
1244                       
1245                        lowestCost = candidateCost;
1246                }
1247        }       
1248
1249#ifdef _DEBUG
1250        Debug << "plane lowest cost: " << lowestCost << endl;
1251#endif
1252        return bestPlane;
1253}
1254
1255int BspTree::GetNextCandidateIdx(int currentIdx, PolygonContainer &polys)
1256{
1257        const int candidateIdx = Random(currentIdx --);
1258
1259        // swap candidates to avoid testing same plane 2 times
1260        std::swap(polys[currentIdx], polys[candidateIdx]);
1261       
1262        return currentIdx;
1263        //return Random((int)polys.size());
1264}
1265
1266float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
1267                                                          const PolygonContainer &polys) const
1268{
1269        float val = 0;
1270
1271        float sumBalancedPolys = 0;
1272        float sumSplits = 0;
1273        float sumPolyArea = 0;
1274        float sumBalancedViewCells = 0;
1275        float sumBlockedRays = 0;
1276        float totalBlockedRays = 0;
1277        //float totalArea = 0;
1278        int totalViewCells = 0;
1279
1280        // need three unique ids for each type of view cell
1281        // for balanced view cells criterium
1282        ViewCell::NewMail();
1283        const int backId = ViewCell::sMailId;
1284        ViewCell::NewMail();
1285        const int frontId = ViewCell::sMailId;
1286        ViewCell::NewMail();
1287        const int frontAndBackId = ViewCell::sMailId;
1288
1289        PolygonContainer::const_iterator it, it_end = polys.end();
1290
1291        for (it = polys.begin(); it != it_end; ++ it)
1292        {
1293                const int classification = (*it)->ClassifyPlane(candidatePlane);
1294
1295                if (mSplitPlaneStrategy & BALANCED_POLYS)
1296                        sumBalancedPolys += sBalancedPolysTable[classification];
1297               
1298                if (mSplitPlaneStrategy & LEAST_SPLITS)
1299                        sumSplits += sLeastPolySplitsTable[classification];
1300
1301                if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
1302                {
1303                        if (classification == Polygon3::COINCIDENT)
1304                                sumPolyArea += (*it)->GetArea();
1305                        //totalArea += area;
1306                }
1307               
1308                if (mSplitPlaneStrategy & BLOCKED_RAYS)
1309                {
1310                        const float blockedRays = (float)(*it)->mPiercingRays.size();
1311               
1312                        if (classification == Polygon3::COINCIDENT)
1313                                sumBlockedRays += blockedRays;
1314                       
1315                        totalBlockedRays += blockedRays;
1316                }
1317
1318                // assign view cells to back or front according to classificaion
1319                if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
1320                {
1321                        MeshInstance *viewCell = (*it)->mParent;
1322               
1323                        // assure that we only count a view cell
1324                        // once for the front and once for the back side of the plane
1325                        if (classification == Polygon3::FRONT_SIDE)
1326                        {
1327                                if ((viewCell->mMailbox != frontId) &&
1328                                        (viewCell->mMailbox != frontAndBackId))
1329                                {
1330                                        sumBalancedViewCells += 1.0;
1331
1332                                        if (viewCell->mMailbox != backId)
1333                                                viewCell->mMailbox = frontId;
1334                                        else
1335                                                viewCell->mMailbox = frontAndBackId;
1336                                       
1337                                        ++ totalViewCells;
1338                                }
1339                        }
1340                        else if (classification == Polygon3::BACK_SIDE)
1341                        {
1342                                if ((viewCell->mMailbox != backId) &&
1343                                    (viewCell->mMailbox != frontAndBackId))
1344                                {
1345                                        sumBalancedViewCells -= 1.0;
1346
1347                                        if (viewCell->mMailbox != frontId)
1348                                                viewCell->mMailbox = backId;
1349                                        else
1350                                                viewCell->mMailbox = frontAndBackId;
1351
1352                                        ++ totalViewCells;
1353                                }
1354                        }
1355                }
1356        }
1357
1358        const float polysSize = (float)polys.size() + Limits::Small;
1359
1360        // all values should be approx. between 0 and 1 so they can be combined
1361        // and scaled with the factors according to their importance
1362        if (mSplitPlaneStrategy & BALANCED_POLYS)
1363                val += mBalancedPolysFactor * fabs(sumBalancedPolys) / polysSize;
1364       
1365        if (mSplitPlaneStrategy & LEAST_SPLITS) 
1366                val += mLeastSplitsFactor * sumSplits / polysSize;
1367
1368        if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
1369                // HACK: polys.size should be total area so scaling is between 0 and 1
1370                val += mLargestPolyAreaFactor * (float)polys.size() / sumPolyArea;
1371
1372        if (mSplitPlaneStrategy & BLOCKED_RAYS)
1373                if (totalBlockedRays != 0)
1374                        val += mBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;
1375
1376        if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
1377                val += mBalancedViewCellsFactor * fabs(sumBalancedViewCells) /
1378                        ((float)totalViewCells + Limits::Small);
1379       
1380        return val;
1381}
1382
1383bool BspTree::BoundRay(const Ray &ray, float &minT, float &maxT) const
1384{
1385        maxT = 1e6;
1386        minT = 0;
1387       
1388        // test with tree bounding box
1389        if (!mBox.GetMinMaxT(ray, &minT, &maxT))
1390                return false;
1391
1392        if (minT < 0) // start ray from origin
1393                minT = 0;
1394
1395        // bound ray or line segment
1396        if ((ray.GetType() == Ray::LOCAL_RAY) &&
1397            !ray.intersections.empty() &&
1398                (ray.intersections[0].mT <= maxT))
1399        {
1400                maxT = ray.intersections[0].mT;
1401        }
1402
1403        return true;
1404}
1405
1406float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
1407                                                          const BoundedRayContainer &rays,
1408                                                          const int pvs,
1409                                                          const float area,
1410                                                          const BspNodeGeometry &cell) const
1411{
1412        float val = 0;
1413
1414        float sumBalancedRays = 0;
1415        float sumRaySplits = 0;
1416       
1417        int backId = 0;
1418        int frontId = 0;
1419        int frontAndBackId = 0;
1420
1421        int frontPvs = 0;
1422        int backPvs = 0;
1423
1424        // probability that view point lies in child
1425        float pOverall = 0;
1426        float pFront = 0;
1427        float pBack = 0;
1428
1429        if (mSplitPlaneStrategy & PVS)
1430        {
1431                // create three unique ids for pvs heuristics
1432                Intersectable::NewMail(); backId = ViewCell::sMailId;
1433                Intersectable::NewMail(); frontId = ViewCell::sMailId;
1434                Intersectable::NewMail(); frontAndBackId = ViewCell::sMailId;
1435
1436                if (mPvsUseArea) // use front and back cell areas to approximate volume
1437                {       
1438                        // construct child geometry with regard to the candidate split plane
1439                        BspNodeGeometry frontCell;
1440                        BspNodeGeometry backCell;
1441
1442                        cell.SplitGeometry(frontCell, backCell, *this, candidatePlane);
1443               
1444                        pFront = frontCell.GetArea();
1445                        pBack = backCell.GetArea();
1446
1447                        pOverall = area;
1448                }
1449        }
1450                       
1451        BoundedRayContainer::const_iterator rit, rit_end = rays.end();
1452
1453        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1454        {
1455                Ray *ray = (*rit)->mRay;
1456                const float minT = (*rit)->mMinT;
1457                const float maxT = (*rit)->mMaxT;
1458
1459                Vector3 entP, extP;
1460
1461                const int cf =
1462                        ray->ClassifyPlane(candidatePlane, minT, maxT, entP, extP);
1463
1464                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1465                {
1466                        sumBalancedRays += sBalancedRaysTable[cf];
1467                }
1468               
1469                if (mSplitPlaneStrategy & BALANCED_RAYS)
1470                {
1471                        sumRaySplits += sLeastRaySplitsTable[cf];
1472                }
1473
1474                if (mSplitPlaneStrategy & PVS)
1475                {
1476                        if (!ray->intersections.empty())
1477                        {
1478                                // in case the ray intersects an objcrs
1479                                // assure that we only count a object
1480                                // once for the front and once for the back side of the plane
1481                                IncPvs(*ray->intersections[0].mObject, frontPvs, backPvs,
1482                                           cf, frontId, backId, frontAndBackId);
1483                        }
1484
1485                        // the source object in the origin of the ray
1486                        if (ray->sourceObject.mObject)
1487                        {
1488                                IncPvs(*ray->sourceObject.mObject, frontPvs, backPvs,
1489                                           cf, frontId, backId, frontAndBackId);
1490                        }
1491
1492                        if (!mPvsUseArea) // use front and back cell areas to approximate volume
1493                        {
1494                                float len = Distance(entP, extP);
1495                                pOverall += len;
1496
1497                                // use length of rays to approximate volume
1498                                switch (cf)
1499                                {
1500                                        case Ray::COINCIDENT:
1501                                                pBack += len;
1502                                                pFront += len;                                         
1503                                                break;
1504                                        case Ray::BACK:
1505                                                pBack += len;
1506                                                break;
1507                                        case Ray::FRONT:
1508                                                pFront += len;
1509                                                break;
1510                                        case Ray::FRONT_BACK:
1511                                                {
1512                                                        // find intersection of ray segment with plane
1513                                                        const Vector3 extp = ray->Extrap(maxT);
1514                                                        const float t = candidatePlane.FindT(ray->GetLoc(), extp);
1515                                       
1516                                                        const float newT = t * maxT;
1517                                                        float newLen = Distance(ray->Extrap(newT), extp);
1518
1519                                                        pFront += len - newLen;
1520                                                        pBack += newLen;
1521                                                }
1522                                                break;
1523                                        case Ray::BACK_FRONT:
1524                                                {
1525                                                        // find intersection of ray segment with plane
1526                                                        const Vector3 extp = ray->Extrap(maxT);
1527                                                        const float t = candidatePlane.FindT(ray->GetLoc(), extp);
1528                                       
1529                                                        const float newT = t * maxT;
1530                                                        float newLen = Distance(ray->Extrap(newT), extp);
1531
1532                                                        pFront += len;
1533                                                        pBack += len - newLen;
1534                                                }
1535                                                break;
1536                                        default:
1537                                                Debug << "Should not come here" << endl;
1538                                                break;
1539                                }
1540                        }
1541                }
1542        }
1543
1544        const float raysSize = (float)rays.size() + Limits::Small;
1545        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1546                val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
1547
1548        if (mSplitPlaneStrategy & BALANCED_RAYS)
1549                val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
1550
1551        float denom = pOverall * (float)pvs * 2.0f + Limits::Small;
1552        if ((mSplitPlaneStrategy & PVS) && area && pvs)
1553        {
1554                val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
1555
1556                // give penalty to unbalanced split
1557                if (0)
1558                if (((pFront * 0.2 + Limits::Small) > pBack) || (pFront < (pBack * 0.2 + Limits::Small)))
1559                        val += 0.5;
1560        }
1561
1562#ifdef _DEBUG
1563        Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
1564                  << " frontpvs: " << frontPvs << " pFront: " << pFront
1565                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
1566#endif
1567        return val;
1568}
1569
1570void BspTree::IncPvs(Intersectable &obj,
1571                                         int &frontPvs,
1572                                         int &backPvs,
1573                                         const int cf,
1574                                         const int frontId,
1575                                         const int backId,
1576                                         const int frontAndBackId) const
1577{
1578        // TODO: does this really belong to no pvs?
1579        //if (cf == Ray::COINCIDENT) return;
1580
1581        if (cf == Ray::FRONT)
1582        {
1583                if ((obj.mMailbox != frontId) &&
1584                        (obj.mMailbox != frontAndBackId))
1585                {
1586                        ++ frontPvs;
1587
1588                        if (obj.mMailbox != backId)
1589                                obj.mMailbox = frontId;
1590                        else
1591                                obj.mMailbox = frontAndBackId;                                 
1592                }
1593        }
1594        else if (cf == Ray::BACK)
1595        {
1596                if ((obj.mMailbox != backId) &&
1597                        (obj.mMailbox != frontAndBackId))
1598                {
1599                        ++ backPvs;
1600
1601                        if (obj.mMailbox != frontId)
1602                                obj.mMailbox = backId;
1603                        else
1604                                obj.mMailbox = frontAndBackId;
1605                }
1606        }
1607        // object belongs to both PVS
1608        else if ((cf == Ray::FRONT_BACK) || (cf == Ray::BACK_FRONT) ||(cf == Ray::COINCIDENT))
1609        {
1610                if (obj.mMailbox !=  frontAndBackId)
1611                {
1612                        if (obj.mMailbox != frontId)
1613                                ++ frontPvs;
1614                        if (obj.mMailbox != backId)
1615                                ++ backPvs;
1616               
1617                        obj.mMailbox = frontAndBackId;
1618                }
1619        }
1620}
1621
1622float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
1623                                                          BspTraversalData &data) const
1624{
1625        float val = 0;
1626
1627        if (mSplitPlaneStrategy & VERTICAL_AXIS)
1628        {
1629                Vector3 tinyAxis(0,0,0); tinyAxis[mBox.Size().TinyAxis()] = 1.0f;
1630                // we put a penalty on the dot product between the "tiny" vertical axis
1631                // and the split plane axis
1632                val += mVerticalSplitsFactor *
1633                           fabs(DotProd(candidatePlane.mNormal, tinyAxis));
1634        }
1635
1636        // the following criteria loop over all polygons to find the cost value
1637        if ((mSplitPlaneStrategy & BALANCED_POLYS)      ||
1638                (mSplitPlaneStrategy & LEAST_SPLITS)        ||
1639                (mSplitPlaneStrategy & LARGEST_POLY_AREA)   ||
1640                (mSplitPlaneStrategy & BALANCED_VIEW_CELLS) ||
1641                (mSplitPlaneStrategy & BLOCKED_RAYS))
1642        {
1643                val += SplitPlaneCost(candidatePlane, *data.mPolygons);
1644        }
1645
1646        // the following criteria loop over all rays to find the cost value
1647        if ((mSplitPlaneStrategy & BALANCED_RAYS)      ||
1648                (mSplitPlaneStrategy & LEAST_RAY_SPLITS)   ||
1649                (mSplitPlaneStrategy & PVS))
1650        {
1651                val += SplitPlaneCost(candidatePlane, *data.mRays, data.mPvs,
1652                                                          data.mArea, *data.mGeometry);
1653        }
1654
1655        // return linear combination of the sums
1656        return val;
1657}
1658
1659void BspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
1660{
1661        stack<BspNode *> nodeStack;
1662        nodeStack.push(mRoot);
1663 
1664        while (!nodeStack.empty())
1665        {
1666                BspNode *node = nodeStack.top();
1667   
1668                nodeStack.pop();
1669   
1670                if (node->IsLeaf())
1671                {
1672                        BspLeaf *leaf = (BspLeaf *)node;               
1673                        leaves.push_back(leaf);
1674                }
1675                else
1676                {
1677                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1678
1679                        nodeStack.push(interior->GetBack());
1680                        nodeStack.push(interior->GetFront());
1681                }
1682        }
1683}
1684
1685AxisAlignedBox3 BspTree::GetBoundingBox() const
1686{
1687        return mBox;
1688}
1689
1690BspNode *BspTree::GetRoot() const
1691{
1692        return mRoot;
1693}
1694
1695void BspTree::EvaluateLeafStats(const BspTraversalData &data)
1696{
1697        // the node became a leaf -> evaluate stats for leafs
1698        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
1699
1700        if (data.mDepth >= mTermMaxDepth)
1701                ++ mStat.maxDepthNodes;
1702       
1703        // store maximal and minimal depth
1704        if (data.mDepth > mStat.maxDepth)
1705                mStat.maxDepth = data.mDepth;
1706
1707        if (data.mDepth < mStat.minDepth)
1708                mStat.minDepth = data.mDepth;
1709
1710        // accumulate depth to compute average depth
1711        mStat.accumDepth += data.mDepth;
1712       
1713#ifdef _DEBUG
1714        Debug << "BSP stats: "
1715                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1716                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1717                  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
1718                  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << mTermMaxPolygons << "), "
1719                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMaxRays << "), "
1720                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
1721                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1722#endif
1723}
1724
1725int BspTree::CastRay(Ray &ray)
1726{
1727        int hits = 0;
1728 
1729        stack<BspRayTraversalData> tStack;
1730 
1731        float maxt, mint;
1732
1733        if (!BoundRay(ray, mint, maxt))
1734                return 0;
1735
1736        Intersectable::NewMail();
1737
1738        Vector3 entp = ray.Extrap(mint);
1739        Vector3 extp = ray.Extrap(maxt);
1740 
1741        BspNode *node = mRoot;
1742        BspNode *farChild = NULL;
1743       
1744        while (1)
1745        {
1746                if (!node->IsLeaf())
1747                {
1748                        BspInterior *in = (BspInterior *) node;
1749                       
1750                        Plane3 *splitPlane = in->GetPlane();
1751
1752                        int entSide = splitPlane->Side(entp);
1753                        int extSide = splitPlane->Side(extp);
1754
1755                        Vector3 intersection;
1756
1757                        if (entSide < 0)
1758                        {
1759                                node = in->GetBack();
1760
1761                                if(extSide <= 0) // plane does not split ray => no far child
1762                                        continue;
1763                                       
1764                                farChild = in->GetFront(); // plane splits ray
1765
1766                        } else if (entSide > 0)
1767                        {
1768                                node = in->GetFront();
1769
1770                                if (extSide >= 0) // plane does not split ray => no far child
1771                                        continue;
1772
1773                                farChild = in->GetBack(); // plane splits ray                   
1774                        }
1775                        else // ray and plane are coincident
1776                        {
1777                                // WHAT TO DO IN THIS CASE ?
1778                                //break;
1779                                node = in->GetFront();
1780                                continue;
1781                        }
1782
1783                        // push data for far child
1784                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
1785
1786                        // find intersection of ray segment with plane
1787                        float t;
1788                        extp = splitPlane->FindIntersection(ray.GetLoc(), extp, &t);
1789                        maxt *= t;
1790                       
1791                } else // reached leaf => intersection with view cell
1792                {
1793                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1794     
1795                        if (!leaf->mViewCell->Mailed())
1796                        {
1797                                ray.bspIntersections.push_back(Ray::BspIntersection(maxt, leaf));
1798                                leaf->mViewCell->Mail();
1799                                ++ hits;
1800                        }
1801                       
1802                        //-- fetch the next far child from the stack
1803                        if (tStack.empty())
1804                                break;
1805     
1806                        entp = extp;
1807                        mint = maxt; // NOTE: need this?
1808
1809                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
1810                                break;
1811
1812                        BspRayTraversalData &s = tStack.top();
1813
1814                        node = s.mNode;
1815                        extp = s.mExitPoint;
1816                        maxt = s.mMaxT;
1817
1818                        tStack.pop();
1819                }
1820        }
1821
1822        return hits;
1823}
1824
1825bool BspTree::Export(const string filename)
1826{
1827        Exporter *exporter = Exporter::GetExporter(filename);
1828
1829        if (exporter)
1830        {
1831                exporter->ExportBspTree(*this);
1832                return true;
1833        }       
1834
1835        return false;
1836}
1837
1838void BspTree::CollectViewCells(ViewCellContainer &viewCells) const
1839{
1840        stack<BspNode *> nodeStack;
1841        nodeStack.push(mRoot);
1842
1843        ViewCell::NewMail();
1844
1845        while (!nodeStack.empty())
1846        {
1847                BspNode *node = nodeStack.top();
1848                nodeStack.pop();
1849
1850                if (node->IsLeaf())
1851                {
1852                        ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->mViewCell;
1853
1854                        if (!viewCell->Mailed())
1855                        {
1856                                viewCell->Mail();
1857                                viewCells.push_back(viewCell);
1858                        }
1859                }
1860                else
1861                {
1862                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1863
1864                        nodeStack.push(interior->mFront);
1865                        nodeStack.push(interior->mBack);
1866                }
1867        }
1868}
1869
1870void BspTree::EvaluateViewCellsStats(BspViewCellsStatistics &stat) const
1871{
1872        stat.Reset();
1873
1874        stack<BspNode *> nodeStack;
1875        nodeStack.push(mRoot);
1876
1877        ViewCell::NewMail();
1878
1879        // exclude root cell
1880        mRootCell->Mail();
1881
1882        while (!nodeStack.empty())
1883        {
1884                BspNode *node = nodeStack.top();
1885                nodeStack.pop();
1886
1887                if (node->IsLeaf())
1888                {
1889                        ++ stat.bspLeaves;
1890
1891                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->mViewCell;
1892
1893                        if (!viewCell->Mailed())
1894                        {
1895                                viewCell->Mail();
1896                               
1897                                ++ stat.viewCells;
1898                                const int pvsSize = viewCell->GetPvs().GetSize();
1899
1900                stat.pvs += pvsSize;
1901
1902                                if (pvsSize < 1)
1903                                        ++ stat.emptyPvs;
1904
1905                                if (pvsSize > stat.maxPvs)
1906                                        stat.maxPvs = pvsSize;
1907
1908                                if (pvsSize < stat.minPvs)
1909                                        stat.minPvs = pvsSize;
1910
1911                                if ((int)viewCell->mBspLeaves.size() > stat.maxBspLeaves)
1912                                        stat.maxBspLeaves = (int)viewCell->mBspLeaves.size();           
1913                        }
1914                }
1915                else
1916                {
1917                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1918
1919                        nodeStack.push(interior->mFront);
1920                        nodeStack.push(interior->mBack);
1921                }
1922        }
1923}
1924
1925bool BspTree::MergeViewCells(BspLeaf *front, BspLeaf *back) const
1926{
1927        BspViewCell *viewCell =
1928                dynamic_cast<BspViewCell *>(ViewCell::Merge(*front->mViewCell, *back->mViewCell));
1929       
1930        if (!viewCell)
1931                return false;
1932
1933        // change view cells of all leaves associated with the
1934        // previous view cells
1935
1936        BspViewCell *fVc = front->mViewCell;
1937        BspViewCell *bVc = back->mViewCell;
1938
1939        vector<BspLeaf *> fLeaves = fVc->mBspLeaves;
1940        vector<BspLeaf *> bLeaves = bVc->mBspLeaves;
1941
1942        vector<BspLeaf *>::const_iterator it;
1943       
1944        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
1945        {
1946                (*it)->SetViewCell(viewCell);
1947                viewCell->mBspLeaves.push_back(*it);
1948        }
1949        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
1950        {
1951                (*it)->SetViewCell(viewCell);
1952                viewCell->mBspLeaves.push_back(*it);
1953        }
1954       
1955        DEL_PTR(fVc);
1956        DEL_PTR(bVc);
1957
1958        return true;
1959}
1960
1961bool BspTree::ShouldMerge(BspLeaf *front, BspLeaf *back) const
1962{
1963        ViewCell *fvc = front->mViewCell;
1964        ViewCell *bvc = back->mViewCell;
1965
1966        if ((fvc == mRootCell) || (bvc == mRootCell) || (fvc == bvc))
1967                return false;
1968
1969        int fdiff = fvc->GetPvs().Diff(bvc->GetPvs());
1970
1971        if (fvc->GetPvs().GetSize() + fdiff < mMaxPvs)
1972        {
1973                if ((fvc->GetPvs().GetSize() < mMinPvs) ||     
1974                        (bvc->GetPvs().GetSize() < mMinPvs) ||
1975                        ((fdiff < mMinPvsDif) && (bvc->GetPvs().Diff(fvc->GetPvs()) < mMinPvsDif)))
1976                {
1977                        return true;
1978                }
1979        }
1980       
1981        return false;
1982}
1983
1984void BspTree::SetGenerateViewCells(int generateViewCells)
1985{
1986        mGenerateViewCells = generateViewCells;
1987}
1988
1989BspTreeStatistics &BspTree::GetStat()
1990{
1991        return mStat;
1992}
1993
1994float BspTree::AccumulatedRayLength(BoundedRayContainer &rays) const
1995{
1996        float len = 0;
1997
1998        BoundedRayContainer::const_iterator it, it_end = rays.end();
1999
2000        for (it = rays.begin(); it != it_end; ++ it)
2001        {
2002                len += SqrDistance((*it)->mRay->Extrap((*it)->mMinT),
2003                                                   (*it)->mRay->Extrap((*it)->mMaxT));
2004        }
2005
2006        return len;
2007}
2008
2009int BspTree::SplitRays(const Plane3 &plane,
2010                                           BoundedRayContainer &rays,
2011                                           BoundedRayContainer &frontRays,
2012                                           BoundedRayContainer &backRays)
2013{
2014        int splits = 0;
2015
2016        while (!rays.empty())
2017        {
2018                BoundedRay *bRay = rays.back();
2019                Ray *ray = bRay->mRay;
2020                float minT = bRay->mMinT;
2021                float maxT = bRay->mMaxT;
2022
2023                rays.pop_back();
2024       
2025                Vector3 entP, extP;
2026
2027                const int cf =
2028                        ray->ClassifyPlane(plane, minT, maxT, entP, extP);
2029               
2030                // set id to ray classification
2031                ray->SetId(cf);
2032
2033                switch (cf)
2034                {
2035                case Ray::COINCIDENT: // TODO: should really discard ray?
2036                        frontRays.push_back(bRay);
2037                        //DEL_PTR(bRay);
2038                        break;
2039                case Ray::BACK:
2040                        backRays.push_back(bRay);
2041                        break;
2042                case Ray::FRONT:
2043                        frontRays.push_back(bRay);
2044                        break;
2045                case Ray::FRONT_BACK:
2046                        {
2047                                // find intersection of ray segment with plane
2048                                const float t = plane.FindT(ray->GetLoc(), extP);
2049                               
2050                                const float newT = t * maxT;
2051
2052                                frontRays.push_back(new BoundedRay(ray, minT, newT));
2053                                backRays.push_back(new BoundedRay(ray, newT, maxT));
2054
2055                                DEL_PTR(bRay);
2056                        }
2057                        break;
2058                case Ray::BACK_FRONT:
2059                        {
2060                                // find intersection of ray segment with plane
2061                                const float t = plane.FindT(ray->GetLoc(), extP);
2062                                const float newT = t * bRay->mMaxT;
2063
2064                                backRays.push_back(new BoundedRay(ray, bRay->mMinT, newT));
2065                                frontRays.push_back(new BoundedRay(ray, newT, bRay->mMaxT));
2066                                DEL_PTR(bRay);
2067
2068                                ++ splits;
2069                        }
2070                        break;
2071                default:
2072                        Debug << "Should not come here" << endl;
2073                        break;
2074                }
2075        }
2076
2077        return splits;
2078}
2079
2080void BspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
2081{
2082        BspNode *lastNode;
2083        do
2084        {
2085                lastNode = n;
2086
2087                // want to get planes defining geometry of this node => don't take
2088                // split plane of node itself
2089                n = n->GetParent();
2090               
2091                if (n)
2092                {
2093                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
2094                        Plane3 halfSpace = *dynamic_cast<BspInterior *>(interior)->GetPlane();
2095
2096            if (interior->mFront != lastNode)
2097                                halfSpace.ReverseOrientation();
2098
2099                        halfSpaces.push_back(halfSpace);
2100                }
2101        }
2102        while (n);
2103}
2104
2105void BspTree::ConstructGeometry(BspNode *n, BspNodeGeometry &cell) const
2106{
2107        PolygonContainer polys;
2108        ConstructGeometry(n, polys);
2109        cell.mPolys = polys;
2110}
2111
2112void BspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &cell) const
2113{
2114        vector<BspLeaf *> leaves = vc->mBspLeaves;
2115
2116        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
2117
2118        for (it = leaves.begin(); it != it_end; ++ it)
2119                ConstructGeometry(*it, cell);
2120}
2121
2122
2123void BspTree::ConstructGeometry(BspNode *n, PolygonContainer &cell) const
2124{
2125        vector<Plane3> halfSpaces;
2126        ExtractHalfSpaces(n, halfSpaces);
2127
2128        PolygonContainer candidatePolys;
2129
2130        // bounded planes are added to the polygons (reverse polygons
2131        // as they have to be outfacing
2132        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
2133        {
2134                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
2135               
2136                if (p->Valid())
2137                {
2138                        candidatePolys.push_back(p->CreateReversePolygon());
2139                        DEL_PTR(p);
2140                }
2141        }
2142
2143        // add faces of bounding box (also could be faces of the cell)
2144        for (int i = 0; i < 6; ++ i)
2145        {
2146                VertexContainer vertices;
2147       
2148                for (int j = 0; j < 4; ++ j)
2149                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
2150
2151                candidatePolys.push_back(new Polygon3(vertices));
2152        }
2153
2154        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
2155        {
2156                // polygon is split by all other planes
2157                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
2158                {
2159                        if (i == j) // polygon and plane are coincident
2160                                continue;
2161
2162                        VertexContainer splitPts;
2163                        Polygon3 *frontPoly, *backPoly;
2164
2165                        const int cf = candidatePolys[i]->ClassifyPlane(halfSpaces[j]);
2166                       
2167                        switch (cf)
2168                        {
2169                                case Polygon3::SPLIT:
2170                                        frontPoly = new Polygon3();
2171                                        backPoly = new Polygon3();
2172
2173                                        candidatePolys[i]->Split(halfSpaces[j], *frontPoly,
2174                                                                                         *backPoly, splitPts);
2175
2176                                        DEL_PTR(candidatePolys[i]);
2177
2178                                        if (frontPoly->Valid())
2179                                                candidatePolys[i] = frontPoly;
2180                                        else
2181                                                DEL_PTR(frontPoly);
2182
2183                                        DEL_PTR(backPoly);
2184                                        break;
2185                                case Polygon3::BACK_SIDE:
2186                                        DEL_PTR(candidatePolys[i]);
2187                                        break;
2188                                // just take polygon as it is
2189                                case Polygon3::FRONT_SIDE:
2190                                case Polygon3::COINCIDENT:
2191                                default:
2192                                        break;
2193                        }
2194                }
2195               
2196                if (candidatePolys[i])
2197                        cell.push_back(candidatePolys[i]);
2198        }
2199}
2200
2201
2202int BspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
2203                                                   const bool onlyUnmailed) const
2204{
2205        PolygonContainer cell;
2206
2207        ConstructGeometry(n, cell);
2208
2209        stack<BspNode *> nodeStack;
2210        nodeStack.push(mRoot);
2211               
2212        // planes needed to verify that we found neighbor leaf.
2213        vector<Plane3> halfSpaces;
2214        ExtractHalfSpaces(n, halfSpaces);
2215
2216        while (!nodeStack.empty())
2217        {
2218                BspNode *node = nodeStack.top();
2219                nodeStack.pop();
2220
2221                if (node->IsLeaf())
2222                {
2223            if (node != n && (!onlyUnmailed || !node->Mailed()))
2224                        {
2225                                // test all planes of current node if candidate really
2226                                // is neighbour
2227                                PolygonContainer neighborCandidate;
2228                                ConstructGeometry(node, neighborCandidate);
2229                               
2230                                bool isAdjacent = true;
2231                                for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
2232                                {
2233                                        const int cf =
2234                                                Polygon3::ClassifyPlane(neighborCandidate, halfSpaces[i]);
2235
2236                                        if (cf == Polygon3::BACK_SIDE)
2237                                                isAdjacent = false;
2238                                }
2239
2240                                if (isAdjacent)
2241                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2242
2243                                CLEAR_CONTAINER(neighborCandidate);
2244                        }
2245                }
2246                else
2247                {
2248                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2249       
2250                        const int cf = Polygon3::ClassifyPlane(cell, interior->mPlane);
2251
2252                        if (cf == Polygon3::FRONT_SIDE)
2253                                nodeStack.push(interior->mFront);
2254                        else
2255                                if (cf == Polygon3::BACK_SIDE)
2256                                        nodeStack.push(interior->mBack);
2257                                else
2258                                {
2259                                        // random decision
2260                                        nodeStack.push(interior->mBack);
2261                                        nodeStack.push(interior->mFront);
2262                                }
2263                }
2264        }
2265       
2266        CLEAR_CONTAINER(cell);
2267        return (int)neighbors.size();
2268}
2269
2270BspLeaf *BspTree::GetRandomLeaf(const Plane3 &halfspace)
2271{
2272    stack<BspNode *> nodeStack;
2273        nodeStack.push(mRoot);
2274       
2275        int mask = rand();
2276 
2277        while (!nodeStack.empty())
2278        {
2279                BspNode *node = nodeStack.top();
2280                nodeStack.pop();
2281         
2282                if (node->IsLeaf())
2283                {
2284                        return dynamic_cast<BspLeaf *>(node);
2285                }
2286                else
2287                {
2288                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2289                       
2290                        BspNode *next;
2291       
2292                        PolygonContainer cell;
2293
2294                        // todo: not very efficient: constructs full cell everytime
2295                        ConstructGeometry(interior, cell);
2296
2297                        const int cf = Polygon3::ClassifyPlane(cell, halfspace);
2298
2299                        if (cf == Polygon3::BACK_SIDE)
2300                                next = interior->mFront;
2301                        else
2302                                if (cf == Polygon3::FRONT_SIDE)
2303                                        next = interior->mFront;
2304                        else
2305                        {
2306                                // random decision
2307                                if (mask & 1)
2308                                        next = interior->mBack;
2309                                else
2310                                        next = interior->mFront;
2311                                mask = mask >> 1;
2312                        }
2313
2314                        nodeStack.push(next);
2315                }
2316        }
2317       
2318        return NULL;
2319}
2320
2321BspLeaf *BspTree::GetRandomLeaf(const bool onlyUnmailed)
2322{
2323        stack<BspNode *> nodeStack;
2324       
2325        nodeStack.push(mRoot);
2326
2327        int mask = rand();
2328       
2329        while (!nodeStack.empty())
2330        {
2331                BspNode *node = nodeStack.top();
2332                nodeStack.pop();
2333               
2334                if (node->IsLeaf())
2335                {
2336                        if ( (!onlyUnmailed || !node->Mailed()) )
2337                                return dynamic_cast<BspLeaf *>(node);
2338                }
2339                else
2340                {
2341                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2342
2343                        // random decision
2344                        if (mask & 1)
2345                                nodeStack.push(interior->mBack);
2346                        else
2347                                nodeStack.push(interior->mFront);
2348
2349                        mask = mask >> 1;
2350                }
2351        }
2352       
2353        return NULL;
2354}
2355
2356int BspTree::ComputePvsSize(const BoundedRayContainer &rays) const
2357{
2358        int pvsSize = 0;
2359
2360        BoundedRayContainer::const_iterator rit, rit_end = rays.end();
2361
2362        Intersectable::NewMail();
2363
2364        for (rit = rays.begin(); rit != rays.end(); ++ rit)
2365        {
2366                Ray *ray = (*rit)->mRay;
2367               
2368                if (!ray->intersections.empty())
2369                {
2370                        if (!ray->intersections[0].mObject->Mailed())
2371                        {
2372                                ray->intersections[0].mObject->Mail();
2373                                ++ pvsSize;
2374                        }
2375                }
2376                if (ray->sourceObject.mObject)
2377                {
2378                        if (!ray->sourceObject.mObject->Mailed())
2379                        {
2380                                ray->sourceObject.mObject->Mail();
2381                                ++ pvsSize;
2382                        }
2383                }
2384        }
2385
2386        return pvsSize;
2387}
2388
2389/*************************************************************
2390 *            BspNodeGeometry Implementation                 *
2391 *************************************************************/
2392
2393BspNodeGeometry::~BspNodeGeometry()
2394{
2395        CLEAR_CONTAINER(mPolys);
2396}
2397
2398float BspNodeGeometry::GetArea() const
2399{
2400        return Polygon3::GetArea(mPolys);
2401}
2402
2403void BspNodeGeometry::SplitGeometry(BspNodeGeometry &front,
2404                                                                        BspNodeGeometry &back,
2405                                                                        const BspTree &tree,                                             
2406                                                                        const Plane3 &splitPlane) const
2407{       
2408        // get cross section of new polygon
2409        Polygon3 *planePoly = tree.GetBoundingBox().CrossSection(splitPlane);
2410
2411        planePoly = SplitPolygon(planePoly, tree);
2412
2413        //-- plane poly splits all other cell polygons
2414        for (int i = 0; i < (int)mPolys.size(); ++ i)
2415        {
2416                const int cf = mPolys[i]->ClassifyPlane(splitPlane, 0.00001f);
2417                       
2418                // split new polygon with all previous planes
2419                switch (cf)
2420                {
2421                        case Polygon3::SPLIT:
2422                                {
2423                                        Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);
2424
2425                                        Polygon3 *frontPoly = new Polygon3();
2426                                        Polygon3 *backPoly = new Polygon3();
2427                               
2428                                        VertexContainer splitPts;
2429                                               
2430                                        poly->Split(splitPlane, *frontPoly, *backPoly, splitPts);
2431
2432                                        DEL_PTR(poly);
2433
2434                                        if (frontPoly->Valid())
2435                                                front.mPolys.push_back(frontPoly);
2436                                        if (backPoly->Valid())
2437                                                back.mPolys.push_back(backPoly);
2438                                }
2439                               
2440                                break;
2441                        case Polygon3::BACK_SIDE:
2442                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));                     
2443                                break;
2444                        case Polygon3::FRONT_SIDE:
2445                                front.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));     
2446                                break;
2447                        case Polygon3::COINCIDENT:
2448                                //front.mPolys.push_back(CreateReversePolygon(mPolys[i]));
2449                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));
2450                                break;
2451                        default:
2452                                break;
2453                }
2454        }
2455
2456        //-- finally add the new polygon to the child cells
2457        if (planePoly)
2458        {
2459                // add polygon with normal pointing into positive half space to back cell
2460                back.mPolys.push_back(planePoly);
2461                // add polygon with reverse orientation to front cell
2462                front.mPolys.push_back(planePoly->CreateReversePolygon());
2463        }
2464
2465        //Debug << "returning new geometry " << mPolys.size() << " f: " << front.mPolys.size() << " b: " << back.mPolys.size() << endl;
2466        //Debug << "old area " << GetArea() << " f: " << front.GetArea() << " b: " << back.GetArea() << endl;
2467}
2468
2469Polygon3 *BspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
2470                                                                                const BspTree &tree) const
2471{
2472        // polygon is split by all other planes
2473        for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
2474        {
2475                Plane3 plane = mPolys[i]->GetSupportingPlane();
2476
2477                const int cf =
2478                        planePoly->ClassifyPlane(plane, 0.00001f);
2479                       
2480                // split new polygon with all previous planes
2481                switch (cf)
2482                {
2483                        case Polygon3::SPLIT:
2484                                {
2485                                        VertexContainer splitPts;
2486                               
2487                                        Polygon3 *frontPoly = new Polygon3();
2488                                        Polygon3 *backPoly = new Polygon3();
2489
2490                                        planePoly->Split(plane, *frontPoly, *backPoly, splitPts);
2491                                        DEL_PTR(planePoly);
2492
2493                                        if (backPoly->Valid())
2494                                                planePoly = backPoly;
2495                                        else
2496                                                DEL_PTR(backPoly);
2497                                }
2498                                break;
2499                        case Polygon3::FRONT_SIDE:
2500                                DEL_PTR(planePoly);
2501                break;
2502                        // polygon is taken as it is
2503                        case Polygon3::BACK_SIDE:
2504                        case Polygon3::COINCIDENT:
2505                        default:
2506                                break;
2507                }
2508        }
2509
2510        return planePoly;
2511}
2512
2513void BspViewCellsStatistics::Print(ostream &app) const
2514{
2515        app << "===== BspViewCells statistics ===============\n";
2516
2517        app << setprecision(4);
2518
2519        //app << "#N_CTIME  ( Construction time [s] )\n" << Time() << " \n";
2520
2521        app << "#N_OVERALLPVS ( objects in PVS )\n" << pvs << endl;
2522
2523        app << "#N_PMAXPVS ( largest PVS )\n" << maxPvs << endl;
2524
2525        app << "#N_PMINPVS ( smallest PVS )\n" << minPvs << endl;
2526
2527        app << "#N_PAVGPVS ( average PVS )\n" << AvgPvs() << endl;
2528
2529        app << "#N_PEMPTYPVS ( view cells with PVS smaller 2 )\n" << emptyPvs << endl;
2530
2531        app << "#N_VIEWCELLS ( number of view cells)\n" << viewCells << endl;
2532
2533        app << "#N_AVGBSPLEAVES (average number of BSP leaves per view cell )\n" << AvgBspLeaves() << endl;
2534
2535        app << "#N_MAXBSPLEAVES ( maximal number of BSP leaves per view cell )\n" << maxBspLeaves << endl;
2536       
2537        app << "===== END OF BspViewCells statistics ==========\n";
2538}
Note: See TracBrowser for help on using the repository browser.