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

Revision 437, 63.7 KB checked in by mattausch, 19 years ago (diff)

detected leak in BspTree?
added specialised fast view cell bsp tree called VspBspTree?

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