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

Revision 436, 63.8 KB checked in by mattausch, 19 years ago (diff)

bsptree view cells working in VssPreprocessor?

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