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

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