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

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