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

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