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

Revision 404, 61.2 KB checked in by mattausch, 19 years ago (diff)

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