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

Revision 590, 75.0 KB checked in by mattausch, 19 years ago (diff)

implemented some code for merge history loading

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_PMINPROBABILITYAREALEAVES  ( 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
866bool BspTree::TerminationCriteriaMet(const BspTraversalData &data) const
867{
868        return
869                (((int)data.mPolygons->size() <= mTermMinPolys) ||
870                 ((int)data.mRays->size() <= mTermMinRays) ||
871                 (data.mPvs <= mTermMinPvs) ||
872                 (data.mProbability <= mTermMinProbability) ||
873                 (data.mDepth >= mTermMaxDepth) ||
874                 (mStat.Leaves() >= mMaxViewCells) ||
875                 (data.GetAvgRayContribution() > mTermMaxRayContribution));
876}
877
878BspNode *BspTree::Subdivide(BspTraversalStack &tStack, BspTraversalData &tData)
879{
880        //-- terminate traversal 
881        if (TerminationCriteriaMet(tData))             
882        {
883                BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
884       
885                BspViewCell *viewCell;
886
887                // generate new view cell for each leaf
888                if (mGenerateViewCells)
889                        viewCell = new BspViewCell();
890                else
891                        // add view cell to leaf
892                        viewCell = dynamic_cast<BspViewCell *>(tData.mViewCell);
893               
894                leaf->SetViewCell(viewCell);
895                viewCell->mLeaf = leaf;
896
897                if (mUseAreaForPvs)
898                        viewCell->SetArea(tData.mProbability);
899                else
900                        viewCell->SetVolume(tData.mProbability);
901               
902                //-- add pvs
903                if (viewCell != mOutOfBoundsCell)
904                {
905                        int conSamp = 0, sampCon = 0;
906                        AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
907                       
908                        mStat.contributingSamples += conSamp;
909                        mStat.sampleContributions += sampCon;
910                }
911
912                EvaluateLeafStats(tData);
913               
914                //-- clean up
915               
916                // discard polygons
917                CLEAR_CONTAINER(*tData.mPolygons);
918                // discard rays
919                CLEAR_CONTAINER(*tData.mRays);
920
921                DEL_PTR(tData.mPolygons);
922                DEL_PTR(tData.mRays);
923                DEL_PTR(tData.mGeometry);
924
925                return leaf;
926        }
927
928        //-- continue subdivision
929        PolygonContainer coincident;
930       
931        BspTraversalData tFrontData(NULL, new PolygonContainer(), tData.mDepth + 1, mOutOfBoundsCell,
932                                                                new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
933        BspTraversalData tBackData(NULL, new PolygonContainer(), tData.mDepth + 1, mOutOfBoundsCell,
934                                                           new BoundedRayContainer(), 0, 0, new BspNodeGeometry());
935
936        // create new interior node and two leaf nodes
937        BspInterior *interior =
938                SubdivideNode(tData, tFrontData, tBackData, coincident);
939
940
941        // extract view cells from coincident polygons according to plane normal
942    // only if front or back polygons are empty
943        if (!mGenerateViewCells)
944        {
945                ExtractViewCells(tFrontData,
946                                                 tBackData,
947                                                 coincident,
948                                                 interior->mPlane);                     
949        }
950
951        // don't need coincident polygons anymory
952        CLEAR_CONTAINER(coincident);
953
954        // push the children on the stack
955        tStack.push(tFrontData);
956        tStack.push(tBackData);
957
958        // cleanup
959        DEL_PTR(tData.mNode);
960
961        DEL_PTR(tData.mPolygons);
962        DEL_PTR(tData.mRays);
963        DEL_PTR(tData.mGeometry);               
964       
965        return interior;
966}
967
968void BspTree::ExtractViewCells(BspTraversalData &frontData,
969                                                           BspTraversalData &backData,
970                                                           const PolygonContainer &coincident,
971                                                           const Plane3 &splitPlane) const
972{
973        // if not empty, tree is further subdivided => don't have to find view cell
974        bool foundFront = !frontData.mPolygons->empty();
975        bool foundBack = !frontData.mPolygons->empty();
976
977        PolygonContainer::const_iterator it =
978                coincident.begin(), it_end = coincident.end();
979
980        //-- find first view cells in front and back leafs
981        for (; !(foundFront && foundBack) && (it != it_end); ++ it)
982        {
983                if (DotProd((*it)->GetNormal(), splitPlane.mNormal) > 0)
984                {
985                        backData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
986                        foundBack = true;
987                }
988                else
989                {
990                        frontData.mViewCell = dynamic_cast<ViewCell *>((*it)->mParent);
991                        foundFront = true;
992                }
993        }
994}
995
996BspInterior *BspTree::SubdivideNode(BspTraversalData &tData,
997                                                                        BspTraversalData &frontData,
998                                                                        BspTraversalData &backData,
999                                                                        PolygonContainer &coincident)
1000{
1001        mStat.nodes += 2;
1002       
1003        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
1004
1005        long startTime;
1006        if (0)
1007        {
1008                Debug << "*********************" << endl;
1009                startTime = GetTime();
1010        }
1011       
1012        // select subdivision plane
1013        BspInterior *interior =
1014                new BspInterior(SelectPlane(leaf, tData));
1015       
1016        if (0)
1017        {
1018                Debug << "time used for split plane selection: "
1019                          << TimeDiff(startTime, GetTime()) * 1e-3 << "s" << endl;
1020        }
1021#ifdef _DEBUG
1022        Debug << interior << endl;
1023#endif
1024       
1025
1026        if (0)
1027        {
1028                Debug << "number of rays: " << (int)tData.mRays->size() << endl;
1029                Debug << "number of polys: " << (int)tData.mPolygons->size() << endl;
1030
1031                startTime = GetTime();
1032        }       
1033       
1034        // subdivide rays into front and back rays
1035        SplitRays(interior->mPlane, *tData.mRays, *frontData.mRays, *backData.mRays);
1036       
1037        if (0)
1038        {
1039                Debug << "time used for rays splitting: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
1040                startTime = GetTime();
1041        }
1042
1043        // subdivide polygons with plane
1044        mStat.polySplits += SplitPolygons(interior->GetPlane(),
1045                                                                          *tData.mPolygons,
1046                                                                          *frontData.mPolygons,
1047                                                                          *backData.mPolygons,
1048                                                                          coincident);
1049
1050        if (0)
1051        {
1052                Debug << "time used for polygon splitting: " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
1053        }
1054
1055    // compute pvs
1056        frontData.mPvs = ComputePvsSize(*frontData.mRays);
1057        backData.mPvs = ComputePvsSize(*backData.mRays);
1058
1059        // split geometry and compute area
1060        if (1)
1061        {
1062                tData.mGeometry->SplitGeometry(*frontData.mGeometry,
1063                                                                           *backData.mGeometry,
1064                                                                           interior->mPlane,
1065                                                                           mBox,
1066                                                                           mEpsilon);
1067       
1068               
1069                frontData.mProbability = frontData.mGeometry->GetVolume();
1070                backData.mProbability = backData.mGeometry->GetVolume();
1071        }
1072
1073        // compute accumulated ray length
1074        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
1075        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
1076
1077        //-- create front and back leaf
1078
1079        BspInterior *parent = leaf->GetParent();
1080
1081        // replace a link from node's parent
1082        if (!leaf->IsRoot())
1083        {
1084                parent->ReplaceChildLink(leaf, interior);
1085                interior->SetParent(parent);
1086        }
1087        else // new root
1088        {
1089                mRoot = interior;
1090        }
1091
1092        // and setup child links
1093        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
1094       
1095        frontData.mNode = interior->GetFront();
1096        backData.mNode = interior->GetBack();
1097       
1098        //DEL_PTR(leaf);
1099        return interior;
1100}
1101
1102
1103void BspTree::SortSplitCandidates(const PolygonContainer &polys,
1104                                                                  const int axis,
1105                                                                  vector<SortableEntry> &splitCandidates) const
1106{
1107        splitCandidates.clear();
1108
1109        int requestedSize = 2 * (int)polys.size();
1110        // creates a sorted split candidates array 
1111        splitCandidates.reserve(requestedSize);
1112
1113        PolygonContainer::const_iterator it, it_end = polys.end();
1114
1115        AxisAlignedBox3 box;
1116
1117        // insert all queries
1118        for(it = polys.begin(); it != it_end; ++ it)
1119        {
1120                box.Initialize();
1121                box.Include(*(*it));
1122               
1123                splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
1124                splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
1125        }
1126
1127        stable_sort(splitCandidates.begin(), splitCandidates.end());
1128}
1129
1130
1131float BspTree::BestCostRatio(const PolygonContainer &polys,
1132                                                         const AxisAlignedBox3 &box,
1133                                                         const int axis,
1134                                                         float &position,
1135                                                         int &objectsBack,
1136                                                         int &objectsFront) const
1137{
1138        vector<SortableEntry> splitCandidates;
1139
1140        SortSplitCandidates(polys, axis, splitCandidates);
1141       
1142        // go through the lists, count the number of objects left and right
1143        // and evaluate the following cost funcion:
1144        // C = ct_div_ci  + (ol + or)/queries
1145       
1146        int objectsLeft = 0, objectsRight = (int)polys.size();
1147       
1148        float minBox = box.Min(axis);
1149        float maxBox = box.Max(axis);
1150        float boxArea = box.SurfaceArea();
1151 
1152        float minBand = minBox + mSplitBorder * (maxBox - minBox);
1153        float maxBand = minBox + (1.0f - mSplitBorder) * (maxBox - minBox);
1154       
1155        float minSum = 1e20f;
1156        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
1157
1158        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci)
1159        {
1160                switch ((*ci).type)
1161                {
1162                        case SortableEntry::POLY_MIN:
1163                                ++ objectsLeft;
1164                                break;
1165                        case SortableEntry::POLY_MAX:
1166                            -- objectsRight;
1167                                break;
1168                        default:
1169                                break;
1170                }
1171               
1172                if ((*ci).value > minBand && (*ci).value < maxBand)
1173                {
1174                        AxisAlignedBox3 lbox = box;
1175                        AxisAlignedBox3 rbox = box;
1176                        lbox.SetMax(axis, (*ci).value);
1177                        rbox.SetMin(axis, (*ci).value);
1178
1179                        const float sum = objectsLeft * lbox.SurfaceArea() +
1180                                                          objectsRight * rbox.SurfaceArea();
1181     
1182                        if (sum < minSum)
1183                        {
1184                                minSum = sum;
1185                                position = (*ci).value;
1186
1187                                objectsBack = objectsLeft;
1188                                objectsFront = objectsRight;
1189                        }
1190                }
1191        }
1192 
1193        const float oldCost = (float)polys.size();
1194        const float newCost = mAxisAlignedCtDivCi + minSum / boxArea;
1195        const float ratio = newCost / oldCost;
1196
1197
1198#if 0
1199  Debug << "====================" << endl;
1200  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
1201      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
1202#endif
1203  return ratio;
1204}
1205
1206bool BspTree::SelectAxisAlignedPlane(Plane3 &plane,
1207                                                                         const PolygonContainer &polys) const
1208{
1209        AxisAlignedBox3 box;
1210        box.Initialize();
1211       
1212        // create bounding box of region
1213        Polygon3::IncludeInBox(polys, box);
1214       
1215        int objectsBack = 0, objectsFront = 0;
1216        int axis = 0;
1217        float costRatio = MAX_FLOAT;
1218        Vector3 position;
1219
1220        //-- area subdivision
1221        for (int i = 0; i < 3; ++ i)
1222        {
1223                float p = 0;
1224                float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
1225               
1226                if (r < costRatio)
1227                {
1228                        costRatio = r;
1229                        axis = i;
1230                        position = p;
1231                }
1232        }
1233       
1234        if (costRatio >= mMaxCostRatio)
1235                return false;
1236
1237        Vector3 norm(0,0,0); norm[axis] = 1.0f;
1238        plane = Plane3(norm, position);
1239
1240        return true;
1241}
1242
1243
1244Plane3 BspTree::SelectPlane(BspLeaf *leaf, BspTraversalData &data)
1245{
1246        if (data.mPolygons->empty() && data.mRays->empty())
1247        {
1248                Debug << "Warning: No autopartition polygon candidate available\n";
1249       
1250                // return axis aligned split
1251                AxisAlignedBox3 box;
1252                box.Initialize();
1253       
1254                // create bounding box of region
1255                Polygon3::IncludeInBox(*data.mPolygons, box);
1256
1257                const int axis = box.Size().DrivingAxis();
1258                const Vector3 position = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1259
1260                Vector3 norm(0,0,0); norm[axis] = 1.0f;
1261                return Plane3(norm, position);
1262        }
1263       
1264        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
1265                ((int)data.mPolygons->size() > mTermMinPolysForAxisAligned) &&
1266                ((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
1267                ((mTermMinObjectsForAxisAligned < 0) ||
1268                 (Polygon3::ParentObjectsSize(*data.mPolygons) > mTermMinObjectsForAxisAligned)))
1269        {
1270                Plane3 plane;
1271                if (SelectAxisAlignedPlane(plane, *data.mPolygons))
1272                        return plane;
1273        }
1274
1275        // simplest strategy: just take next polygon
1276        if (mSplitPlaneStrategy & RANDOM_POLYGON)
1277        {
1278        if (!data.mPolygons->empty())
1279                {
1280                        Polygon3 *nextPoly =
1281                                (*data.mPolygons)[(int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1))];
1282                        return nextPoly->GetSupportingPlane();
1283                }
1284                else
1285                {
1286                        const int candidateIdx = (int)RandomValue(0, (Real)((int)data.mRays->size() - 1));
1287                        BoundedRay *bRay = (*data.mRays)[candidateIdx];
1288
1289                        Ray *ray = bRay->mRay;
1290                                               
1291                        const Vector3 minPt = ray->Extrap(bRay->mMinT);
1292                        const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
1293
1294                        const Vector3 pt = (maxPt + minPt) * 0.5;
1295
1296                        const Vector3 normal = ray->GetDir();
1297                       
1298                        return Plane3(normal, pt);
1299                }
1300
1301                return Plane3();
1302        }
1303
1304        // use heuristics to find appropriate plane
1305        return SelectPlaneHeuristics(leaf, data);
1306}
1307
1308
1309Plane3 BspTree::ChooseCandidatePlane(const BoundedRayContainer &rays) const
1310{       
1311        const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1312        BoundedRay *bRay = rays[candidateIdx];
1313        Ray *ray = bRay->mRay;
1314
1315        const Vector3 minPt = ray->Extrap(bRay->mMinT);
1316        const Vector3 maxPt = ray->Extrap(bRay->mMaxT);
1317
1318        const Vector3 pt = (maxPt + minPt) * 0.5;
1319
1320        const Vector3 normal = ray->GetDir();
1321                       
1322        return Plane3(normal, pt);
1323}
1324
1325Plane3 BspTree::ChooseCandidatePlane2(const BoundedRayContainer &rays) const
1326{       
1327        Vector3 pt[3];
1328        int idx[3];
1329        int cmaxT = 0;
1330        int cminT = 0;
1331        bool chooseMin = false;
1332
1333        for (int j = 0; j < 3; j ++)
1334        {
1335                idx[j] = (int)RandomValue(0, Real((int)rays.size() * 2 - 1));
1336                               
1337                if (idx[j] >= (int)rays.size())
1338                {
1339                        idx[j] -= (int)rays.size();             
1340                        chooseMin = (cminT < 2);
1341                }
1342                else
1343                        chooseMin = (cmaxT < 2);
1344
1345                BoundedRay *bRay = rays[idx[j]];
1346                pt[j] = chooseMin ? bRay->mRay->Extrap(bRay->mMinT) :
1347                                                        bRay->mRay->Extrap(bRay->mMaxT);
1348        }       
1349                       
1350        return Plane3(pt[0], pt[1], pt[2]);
1351}
1352
1353Plane3 BspTree::ChooseCandidatePlane3(const BoundedRayContainer &rays) const
1354{       
1355        Vector3 pt[3];
1356       
1357        int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1358        int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1359
1360        // check if rays different
1361        if (idx1 == idx2)
1362                idx2 = (idx2 + 1) % (int)rays.size();
1363
1364        const BoundedRay *ray1 = rays[idx1];
1365        const BoundedRay *ray2 = rays[idx2];
1366
1367        // normal vector of the plane parallel to both lines
1368        const Vector3 norm =
1369                Normalize(CrossProd(ray1->mRay->GetDir(), ray2->mRay->GetDir()));
1370
1371        const Vector3 orig1 = ray1->mRay->Extrap(ray1->mMinT);
1372        const Vector3 orig2 = ray2->mRay->Extrap(ray2->mMinT);
1373
1374        // vector from line 1 to line 2
1375        const Vector3 vd = orig1 - orig2;
1376       
1377        // project vector on normal to get distance
1378        const float dist = DotProd(vd, norm);
1379
1380        // point on plane lies halfway between the two planes
1381        const Vector3 planePt = orig1 + norm * dist * 0.5;
1382
1383        return Plane3(norm, planePt);
1384}
1385
1386
1387Plane3 BspTree::SelectPlaneHeuristics(BspLeaf *leaf, BspTraversalData &data)
1388{
1389        float lowestCost = MAX_FLOAT;
1390        Plane3 bestPlane;
1391        // intermediate plane
1392        Plane3 plane;
1393
1394        const int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
1395        int maxIdx = (int)data.mPolygons->size();
1396       
1397        for (int i = 0; i < limit; ++ i)
1398        {
1399                // assure that no index is taken twice
1400                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
1401                               
1402                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
1403
1404                // swap candidate to the end to avoid testing same plane
1405                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
1406                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
1407
1408                // evaluate current candidate
1409                const float candidateCost =
1410                        SplitPlaneCost(poly->GetSupportingPlane(), data);
1411
1412                if (candidateCost < lowestCost)
1413                {
1414                        bestPlane = poly->GetSupportingPlane();
1415                        lowestCost = candidateCost;
1416                }
1417        }
1418       
1419        //-- choose candidate planes extracted from rays
1420        for (int i = 0; i < mMaxRayCandidates; ++ i)
1421        {
1422                plane = ChooseCandidatePlane3(*data.mRays);
1423                const float candidateCost = SplitPlaneCost(plane, data);
1424
1425                if (candidateCost < lowestCost)
1426                {
1427                        bestPlane = plane;     
1428                        lowestCost = candidateCost;
1429                }
1430        }
1431
1432#ifdef _DEBUG
1433        Debug << "plane lowest cost: " << lowestCost << endl;
1434#endif
1435
1436        return bestPlane;
1437}
1438
1439
1440float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
1441                                                          const PolygonContainer &polys) const
1442{
1443        float val = 0;
1444
1445        float sumBalancedPolys = 0;
1446        float sumSplits = 0;
1447        float sumPolyArea = 0;
1448        float sumBalancedViewCells = 0;
1449        float sumBlockedRays = 0;
1450        float totalBlockedRays = 0;
1451        //float totalArea = 0;
1452        int totalViewCells = 0;
1453
1454        // need three unique ids for each type of view cell
1455        // for balanced view cells criterium
1456        ViewCell::NewMail();
1457        const int backId = ViewCell::sMailId;
1458        ViewCell::NewMail();
1459        const int frontId = ViewCell::sMailId;
1460        ViewCell::NewMail();
1461        const int frontAndBackId = ViewCell::sMailId;
1462
1463        bool useRand;;
1464        int limit;
1465
1466        // choose test polyongs randomly if over threshold
1467        if ((int)polys.size() > mMaxTests)
1468        {
1469                useRand = true;
1470                limit = mMaxTests;
1471        }
1472        else
1473        {
1474                useRand = false;
1475                limit = (int)polys.size();
1476        }
1477
1478        for (int i = 0; i < limit; ++ i)
1479        {
1480                const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;
1481
1482                Polygon3 *poly = polys[testIdx];
1483
1484        const int classification =
1485                        poly->ClassifyPlane(candidatePlane, mEpsilon);
1486
1487                if (mSplitPlaneStrategy & BALANCED_POLYS)
1488                        sumBalancedPolys += sBalancedPolysTable[classification];
1489               
1490                if (mSplitPlaneStrategy & LEAST_SPLITS)
1491                        sumSplits += sLeastPolySplitsTable[classification];
1492
1493                if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
1494                {
1495                        if (classification == Polygon3::COINCIDENT)
1496                                sumPolyArea += poly->GetArea();
1497                        //totalArea += area;
1498                }
1499               
1500                if (mSplitPlaneStrategy & BLOCKED_RAYS)
1501                {
1502                        const float blockedRays = (float)poly->mPiercingRays.size();
1503               
1504                        if (classification == Polygon3::COINCIDENT)
1505                                sumBlockedRays += blockedRays;
1506                       
1507                        totalBlockedRays += blockedRays;
1508                }
1509
1510                // assign view cells to back or front according to classificaion
1511                if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
1512                {
1513                        MeshInstance *viewCell = poly->mParent;
1514               
1515                        // assure that we only count a view cell
1516                        // once for the front and once for the back side of the plane
1517                        if (classification == Polygon3::FRONT_SIDE)
1518                        {
1519                                if ((viewCell->mMailbox != frontId) &&
1520                                        (viewCell->mMailbox != frontAndBackId))
1521                                {
1522                                        sumBalancedViewCells += 1.0;
1523
1524                                        if (viewCell->mMailbox != backId)
1525                                                viewCell->mMailbox = frontId;
1526                                        else
1527                                                viewCell->mMailbox = frontAndBackId;
1528                                       
1529                                        ++ totalViewCells;
1530                                }
1531                        }
1532                        else if (classification == Polygon3::BACK_SIDE)
1533                        {
1534                                if ((viewCell->mMailbox != backId) &&
1535                                    (viewCell->mMailbox != frontAndBackId))
1536                                {
1537                                        sumBalancedViewCells -= 1.0;
1538
1539                                        if (viewCell->mMailbox != frontId)
1540                                                viewCell->mMailbox = backId;
1541                                        else
1542                                                viewCell->mMailbox = frontAndBackId;
1543
1544                                        ++ totalViewCells;
1545                                }
1546                        }
1547                }
1548        }
1549
1550        const float polysSize = (float)polys.size() + Limits::Small;
1551
1552        // all values should be approx. between 0 and 1 so they can be combined
1553        // and scaled with the factors according to their importance
1554        if (mSplitPlaneStrategy & BALANCED_POLYS)
1555                val += mBalancedPolysFactor * fabs(sumBalancedPolys) / polysSize;
1556       
1557        if (mSplitPlaneStrategy & LEAST_SPLITS) 
1558                val += mLeastSplitsFactor * sumSplits / polysSize;
1559
1560        if (mSplitPlaneStrategy & LARGEST_POLY_AREA)
1561                // HACK: polys.size should be total area so scaling is between 0 and 1
1562                val += mLargestPolyAreaFactor * (float)polys.size() / sumPolyArea;
1563
1564        if (mSplitPlaneStrategy & BLOCKED_RAYS)
1565                if (totalBlockedRays != 0)
1566                        val += mBlockedRaysFactor * (totalBlockedRays - sumBlockedRays) / totalBlockedRays;
1567
1568        if (mSplitPlaneStrategy & BALANCED_VIEW_CELLS)
1569                val += mBalancedViewCellsFactor * fabs(sumBalancedViewCells) /
1570                        ((float)totalViewCells + Limits::Small);
1571       
1572        return val;
1573}
1574
1575
1576inline void BspTree::GenerateUniqueIdsForPvs()
1577{
1578        Intersectable::NewMail(); sBackId = ViewCell::sMailId;
1579        Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
1580        Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
1581}
1582
1583
1584float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
1585                                                          const BoundedRayContainer &rays,
1586                                                          const int pvs,
1587                                                          const float probability,
1588                                                          const BspNodeGeometry &cell) const
1589{
1590        float val = 0;
1591
1592        float sumBalancedRays = 0;
1593        float sumRaySplits = 0;
1594
1595        int frontPvs = 0;
1596        int backPvs = 0;
1597
1598        // probability that view point lies in child
1599        float pOverall = 0;
1600        float pFront = 0;
1601        float pBack = 0;
1602
1603        const bool pvsUseLen = false;
1604
1605        if (mSplitPlaneStrategy & PVS)
1606        {
1607                // create unique ids for pvs heuristics
1608                GenerateUniqueIdsForPvs();
1609
1610                // construct child geometry with regard to the candidate split plane
1611                BspNodeGeometry frontCell;
1612                BspNodeGeometry backCell;
1613
1614                cell.SplitGeometry(frontCell,
1615                                                   backCell,
1616                                                   candidatePlane,
1617                                                   mBox,
1618                                                   mEpsilon);
1619
1620                if (mUseAreaForPvs)
1621                {
1622                        pFront = frontCell.GetArea();
1623                        pBack = backCell.GetArea();
1624                }
1625                else
1626                {
1627                        pFront = frontCell.GetVolume();
1628                        pBack = backCell.GetVolume();
1629                }
1630               
1631
1632                pOverall = probability;
1633        }
1634                       
1635        bool useRand;
1636        int limit;
1637
1638        // choose test polyongs randomly if over threshold
1639        if ((int)rays.size() > mMaxTests)
1640        {
1641                useRand = true;
1642                limit = mMaxTests;
1643        }
1644        else
1645        {
1646                useRand = false;
1647                limit = (int)rays.size();
1648        }
1649
1650        for (int i = 0; i < limit; ++ i)
1651        {
1652                const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;
1653       
1654                BoundedRay *bRay = rays[testIdx];
1655
1656                Ray *ray = bRay->mRay;
1657                const float minT = bRay->mMinT;
1658                const float maxT = bRay->mMaxT;
1659
1660                Vector3 entP, extP;
1661
1662                const int cf =
1663                        ray->ClassifyPlane(candidatePlane, minT, maxT, entP, extP);
1664
1665                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1666                {
1667                        sumBalancedRays += sBalancedRaysTable[cf];
1668                }
1669               
1670                if (mSplitPlaneStrategy & BALANCED_RAYS)
1671                {
1672                        sumRaySplits += sLeastRaySplitsTable[cf];
1673                }
1674
1675                if (mSplitPlaneStrategy & PVS)
1676                {
1677                        // in case the ray intersects an object
1678                        // assure that we only count the object
1679                        // once for the front and once for the back side of the plane
1680                       
1681                        // add the termination object
1682                        if (!ray->intersections.empty())
1683                                AddObjToPvs(ray->intersections[0].mObject, cf, frontPvs, backPvs);
1684                       
1685                        // add the source object
1686                        AddObjToPvs(ray->sourceObject.mObject, cf, frontPvs, backPvs);
1687                }
1688        }
1689
1690        const float raysSize = (float)rays.size() + Limits::Small;
1691
1692        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1693                val += mLeastRaySplitsFactor * sumRaySplits / raysSize;
1694
1695        if (mSplitPlaneStrategy & BALANCED_RAYS)
1696                val += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
1697
1698        const float denom = pOverall * (float)pvs * 2.0f + Limits::Small;
1699
1700        if (mSplitPlaneStrategy & PVS)
1701        {
1702                val += mPvsFactor * (frontPvs * pFront + (backPvs * pBack)) / denom;
1703
1704                // give penalty to unbalanced split
1705                if (0)
1706                if (((pFront * 0.2 + Limits::Small) > pBack) ||
1707                        (pFront < (pBack * 0.2 + Limits::Small)))
1708                        val += 0.5;
1709        }
1710
1711       
1712#ifdef _DEBUG
1713        Debug << "totalpvs: " << pvs << " ptotal: " << pOverall
1714                  << " frontpvs: " << frontPvs << " pFront: " << pFront
1715                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
1716#endif
1717       
1718        return val;
1719}
1720
1721void BspTree::AddObjToPvs(Intersectable *obj,
1722                                                  const int cf,
1723                                                  int &frontPvs,
1724                                                  int &backPvs) const
1725{
1726        if (!obj)
1727                return;
1728        // TODO: does this really belong to no pvs?
1729        //if (cf == Ray::COINCIDENT) return;
1730
1731        // object belongs to both PVS
1732        const bool bothSides = (cf == Ray::FRONT_BACK) ||
1733                                                   (cf == Ray::BACK_FRONT) ||
1734                                                   (cf == Ray::COINCIDENT);
1735
1736        if ((cf == Ray::FRONT) || bothSides)
1737        {
1738                if ((obj->mMailbox != sFrontId) &&
1739                        (obj->mMailbox != sFrontAndBackId))
1740                {
1741                        ++ frontPvs;
1742
1743                        if (obj->mMailbox == sBackId)
1744                                obj->mMailbox = sFrontAndBackId;       
1745                        else
1746                                obj->mMailbox = sFrontId;                                                               
1747                }
1748        }
1749       
1750        if ((cf == Ray::BACK) || bothSides)
1751        {
1752                if ((obj->mMailbox != sBackId) &&
1753                        (obj->mMailbox != sFrontAndBackId))
1754                {
1755                        ++ backPvs;
1756
1757                        if (obj->mMailbox == sFrontId)
1758                                obj->mMailbox = sFrontAndBackId;
1759                        else
1760                                obj->mMailbox = sBackId;                               
1761                }
1762        }
1763}
1764
1765float BspTree::SplitPlaneCost(const Plane3 &candidatePlane,
1766                                                          BspTraversalData &data) const
1767{
1768        float val = 0;
1769
1770        if (mSplitPlaneStrategy & VERTICAL_AXIS)
1771        {
1772                Vector3 tinyAxis(0,0,0); tinyAxis[mBox.Size().TinyAxis()] = 1.0f;
1773                // we put a penalty on the dot product between the "tiny" vertical axis
1774                // and the split plane axis
1775                val += mVerticalSplitsFactor *
1776                           fabs(DotProd(candidatePlane.mNormal, tinyAxis));
1777        }
1778
1779        // the following criteria loop over all polygons to find the cost value
1780        if ((mSplitPlaneStrategy & BALANCED_POLYS)      ||
1781                (mSplitPlaneStrategy & LEAST_SPLITS)        ||
1782                (mSplitPlaneStrategy & LARGEST_POLY_AREA)   ||
1783                (mSplitPlaneStrategy & BALANCED_VIEW_CELLS) ||
1784                (mSplitPlaneStrategy & BLOCKED_RAYS))
1785        {
1786                val += SplitPlaneCost(candidatePlane, *data.mPolygons);
1787        }
1788
1789        // the following criteria loop over all rays to find the cost value
1790        if ((mSplitPlaneStrategy & BALANCED_RAYS)      ||
1791                (mSplitPlaneStrategy & LEAST_RAY_SPLITS)   ||
1792                (mSplitPlaneStrategy & PVS))
1793        {
1794                val += SplitPlaneCost(candidatePlane, *data.mRays, data.mPvs,
1795                                                          data.mProbability, *data.mGeometry);
1796        }
1797
1798        // return linear combination of the sums
1799        return val;
1800}
1801
1802void BspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
1803{
1804        stack<BspNode *> nodeStack;
1805        nodeStack.push(mRoot);
1806 
1807        while (!nodeStack.empty())
1808        {
1809                BspNode *node = nodeStack.top();
1810   
1811                nodeStack.pop();
1812   
1813                if (node->IsLeaf())
1814                {
1815                        BspLeaf *leaf = (BspLeaf *)node;               
1816                        leaves.push_back(leaf);
1817                }
1818                else
1819                {
1820                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1821
1822                        nodeStack.push(interior->GetBack());
1823                        nodeStack.push(interior->GetFront());
1824                }
1825        }
1826}
1827
1828
1829AxisAlignedBox3 BspTree::GetBoundingBox() const
1830{
1831        return mBox;
1832}
1833
1834
1835void BspTree::EvaluateLeafStats(const BspTraversalData &data)
1836{
1837        // the node became a leaf -> evaluate stats for leafs
1838        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
1839
1840        // store maximal and minimal depth
1841        if (data.mDepth > mStat.maxDepth)
1842                mStat.maxDepth = data.mDepth;
1843
1844        if (data.mDepth < mStat.minDepth)
1845                mStat.minDepth = data.mDepth;
1846
1847        // accumulate depth to compute average depth
1848        mStat.accumDepth += data.mDepth;
1849        // accumulate rays to compute rays /  leaf
1850        mStat.accumRays += (int)data.mRays->size();
1851
1852        if (data.mDepth >= mTermMaxDepth)
1853                ++ mStat.maxDepthNodes;
1854
1855        if (data.mPvs < mTermMinPvs)
1856                ++ mStat.minPvsNodes;
1857
1858        if ((int)data.mRays->size() < mTermMinRays)
1859                ++ mStat.minRaysNodes;
1860
1861        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1862                ++ mStat.maxRayContribNodes;
1863       
1864        if (data.mProbability <= mTermMinProbability)
1865                ++ mStat.minProbabilityNodes;
1866
1867#ifdef _DEBUG
1868        Debug << "BSP stats: "
1869                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1870                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1871                  << "Probability: " << data.mProbability << " (min: " << mTermMinProbability << "), "
1872                  << "#polygons: " << (int)data.mPolygons->size() << " (max: " << mTermMinPolys << "), "
1873                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1874                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
1875                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1876#endif
1877}
1878
1879int
1880BspTree::_CastRay(Ray &ray)
1881{
1882        int hits = 0;
1883 
1884        stack<BspRayTraversalData> tStack;
1885 
1886        float maxt, mint;
1887
1888        if (!mBox.GetRaySegment(ray, mint, maxt))
1889                return 0;
1890
1891        Intersectable::NewMail();
1892
1893        Vector3 entp = ray.Extrap(mint);
1894        Vector3 extp = ray.Extrap(maxt);
1895 
1896        BspNode *node = mRoot;
1897        BspNode *farChild = NULL;
1898       
1899        while (1)
1900        {
1901                if (!node->IsLeaf())
1902                {
1903                        BspInterior *in = dynamic_cast<BspInterior *>(node);
1904                       
1905                        Plane3 splitPlane = in->GetPlane();
1906                        const int entSide = splitPlane.Side(entp);
1907                        const int extSide = splitPlane.Side(extp);
1908
1909                        if (entSide < 0)
1910                        {
1911                                node = in->GetBack();
1912
1913                                if(extSide <= 0) // plane does not split ray => no far child
1914                                        continue;
1915                                       
1916                                farChild = in->GetFront(); // plane splits ray
1917
1918                        } else if (entSide > 0)
1919                        {
1920                                node = in->GetFront();
1921
1922                                if (extSide >= 0) // plane does not split ray => no far child
1923                                        continue;
1924
1925                                farChild = in->GetBack(); // plane splits ray                   
1926                        }
1927                        else // ray and plane are coincident
1928                        {
1929                                // WHAT TO DO IN THIS CASE ?
1930                                //break;
1931                                node = in->GetFront();
1932                                continue;
1933                        }
1934
1935                        // push data for far child
1936                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
1937
1938                        // find intersection of ray segment with plane
1939                        float t;
1940                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
1941                        maxt *= t;
1942                       
1943                } else // reached leaf => intersection with view cell
1944                {
1945                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1946     
1947                        if (!leaf->mViewCell->Mailed())
1948                        {
1949                          //                            ray.bspIntersections.push_back(Ray::BspIntersection(maxt, leaf));
1950                                leaf->mViewCell->Mail();
1951                                ++ hits;
1952                        }
1953                       
1954                        //-- fetch the next far child from the stack
1955                        if (tStack.empty())
1956                                break;
1957     
1958                        entp = extp;
1959                        mint = maxt; // NOTE: need this?
1960
1961                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
1962                                break;
1963
1964                        BspRayTraversalData &s = tStack.top();
1965
1966                        node = s.mNode;
1967                        extp = s.mExitPoint;
1968                        maxt = s.mMaxT;
1969
1970                        tStack.pop();
1971                }
1972        }
1973
1974        return hits;
1975}
1976
1977
1978int BspTree::CastLineSegment(const Vector3 &origin,
1979                                                         const Vector3 &termination,
1980                                                         vector<ViewCell *> &viewcells)
1981{
1982        int hits = 0;
1983        stack<BspRayTraversalData> tStack;
1984
1985        float mint = 0.0f, maxt = 1.0f;
1986
1987        Intersectable::NewMail();
1988
1989        Vector3 entp = origin;
1990        Vector3 extp = termination;
1991
1992        BspNode *node = mRoot;
1993        BspNode *farChild = NULL;
1994
1995        while (1)
1996        {
1997                if (!node->IsLeaf()) 
1998                {
1999                        BspInterior *in = dynamic_cast<BspInterior *>(node);
2000               
2001                        Plane3 splitPlane = in->GetPlane();
2002                       
2003                        const int entSide = splitPlane.Side(entp);
2004                        const int extSide = splitPlane.Side(extp);
2005
2006                        if (entSide < 0)
2007                        {
2008                                node = in->GetBack();
2009               
2010                                if(extSide <= 0) // plane does not split ray => no far child
2011                                        continue;
2012               
2013                                farChild = in->GetFront(); // plane splits ray
2014               
2015                        }
2016                        else if (entSide > 0)
2017                        {
2018                                node = in->GetFront();
2019                       
2020                                if (extSide >= 0) // plane does not split ray => no far child
2021                                        continue;
2022                       
2023                                farChild = in->GetBack(); // plane splits ray                   
2024                        }
2025                        else // ray and plane are coincident
2026                        {
2027                                // WHAT TO DO IN THIS CASE ?
2028                                //break;
2029                                node = in->GetFront();
2030                                continue;
2031                        }
2032               
2033                        // push data for far child
2034                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
2035               
2036                        // find intersection of ray segment with plane
2037                        float t;
2038                        extp = splitPlane.FindIntersection(origin, extp, &t);
2039                        maxt *= t; 
2040                }
2041                else
2042                {
2043                        // reached leaf => intersection with view cell
2044                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
2045               
2046                        if (!leaf->mViewCell->Mailed())
2047                        {
2048                                viewcells.push_back(leaf->mViewCell);
2049                                leaf->mViewCell->Mail();
2050                                hits++;
2051                        }
2052               
2053                        //-- fetch the next far child from the stack
2054                        if (tStack.empty())
2055                                break;
2056           
2057                        entp = extp;
2058                        mint = maxt; // NOTE: need this?
2059               
2060                        BspRayTraversalData &s = tStack.top();
2061               
2062                        node = s.mNode;
2063                        extp = s.mExitPoint;
2064                        maxt = s.mMaxT;
2065               
2066                        tStack.pop();
2067                }
2068        }
2069        return hits;
2070}
2071
2072
2073void BspTree::CollectViewCells(ViewCellContainer &viewCells) const
2074{
2075        stack<BspNode *> nodeStack;
2076        nodeStack.push(mRoot);
2077
2078        ViewCell::NewMail();
2079
2080        while (!nodeStack.empty())
2081        {
2082                BspNode *node = nodeStack.top();
2083                nodeStack.pop();
2084
2085                if (node->IsLeaf())
2086                {
2087                        ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->mViewCell;
2088
2089                        if (!viewCell->Mailed())
2090                        {
2091                                viewCell->Mail();
2092                                viewCells.push_back(viewCell);
2093                        }
2094                }
2095                else
2096                {
2097                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2098
2099                        nodeStack.push(interior->GetFront());
2100                        nodeStack.push(interior->GetBack());
2101                }
2102        }
2103}
2104
2105
2106BspTreeStatistics &BspTree::GetStat()
2107{
2108        return mStat;
2109}
2110
2111
2112float BspTree::AccumulatedRayLength(BoundedRayContainer &rays) const
2113{
2114        float len = 0;
2115
2116        BoundedRayContainer::const_iterator it, it_end = rays.end();
2117
2118        for (it = rays.begin(); it != it_end; ++ it)
2119        {
2120                len += SqrDistance((*it)->mRay->Extrap((*it)->mMinT),
2121                                                   (*it)->mRay->Extrap((*it)->mMaxT));
2122        }
2123
2124        return len;
2125}
2126
2127
2128int BspTree::SplitRays(const Plane3 &plane,
2129                                           BoundedRayContainer &rays,
2130                                           BoundedRayContainer &frontRays,
2131                                           BoundedRayContainer &backRays)
2132{
2133        int splits = 0;
2134       
2135        while (!rays.empty())
2136        {
2137                BoundedRay *bRay = rays.back();
2138                Ray *ray = bRay->mRay;
2139                float minT = bRay->mMinT;
2140                float maxT = bRay->mMaxT;
2141
2142                rays.pop_back();
2143       
2144                Vector3 entP, extP;
2145
2146                const int cf =
2147                        ray->ClassifyPlane(plane, minT, maxT, entP, extP);
2148               
2149                // set id to ray classification
2150                ray->SetId(cf);
2151
2152                switch (cf)
2153                {
2154                case Ray::COINCIDENT: // TODO: should really discard ray?
2155                        frontRays.push_back(bRay);
2156                        //DEL_PTR(bRay);
2157                        break;
2158                case Ray::BACK:
2159                        backRays.push_back(bRay);
2160                        break;
2161                case Ray::FRONT:
2162                        frontRays.push_back(bRay);
2163                        break;
2164                case Ray::FRONT_BACK:
2165                        {
2166                                // find intersection of ray segment with plane
2167                                const float t = plane.FindT(ray->GetLoc(), extP);
2168                               
2169                                const float newT = t * maxT;
2170
2171                                frontRays.push_back(new BoundedRay(ray, minT, newT));
2172                                backRays.push_back(new BoundedRay(ray, newT, maxT));
2173
2174                                DEL_PTR(bRay);
2175                        }
2176                        break;
2177                case Ray::BACK_FRONT:
2178                        {
2179                                // find intersection of ray segment with plane
2180                                const float t = plane.FindT(ray->GetLoc(), extP);
2181                                const float newT = t * bRay->mMaxT;
2182
2183                                backRays.push_back(new BoundedRay(ray, minT, newT));
2184                                frontRays.push_back(new BoundedRay(ray, newT, maxT));
2185
2186                                DEL_PTR(bRay);
2187
2188                                ++ splits;
2189                        }
2190                        break;
2191                default:
2192                        Debug << "Should not come here 4" << endl;
2193                        break;
2194                }
2195        }
2196
2197        return splits;
2198}
2199
2200void BspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
2201{
2202        BspNode *lastNode;
2203        do
2204        {
2205                lastNode = n;
2206
2207                // want to get planes defining geometry of this node => don't take
2208                // split plane of node itself
2209                n = n->GetParent();
2210               
2211                if (n)
2212                {
2213                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
2214                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
2215
2216            if (interior->GetFront() != lastNode)
2217                                halfSpace.ReverseOrientation();
2218
2219                        halfSpaces.push_back(halfSpace);
2220                }
2221        }
2222        while (n);
2223}
2224
2225
2226
2227void BspTree::ConstructGeometry(ViewCell *vc,
2228                                                                BspNodeGeometry &vcGeom) const
2229{
2230        ViewCellContainer leaves;
2231        mViewCellsManager->GetViewCellsTree()->CollectLeaves(vc, leaves);
2232
2233        ViewCellContainer::const_iterator it, it_end = leaves.end();
2234
2235        for (it = leaves.begin(); it != it_end; ++ it)
2236        {
2237                BspLeaf *l = dynamic_cast<BspViewCell *>(*it)->mLeaf;
2238                ConstructGeometry(l, vcGeom);
2239        }
2240}
2241
2242
2243void BspTree::SetViewCellsManager(ViewCellsManager *vcm)
2244{
2245        mViewCellsManager = vcm;
2246}
2247
2248
2249void BspTree::ConstructGeometry(BspNode *n, BspNodeGeometry &geom) const
2250{
2251        vector<Plane3> halfSpaces;
2252        ExtractHalfSpaces(n, halfSpaces);
2253
2254        PolygonContainer candidates;
2255
2256        // bounded planes are added to the polygons (reverse polygons
2257        // as they have to be outfacing
2258        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
2259        {
2260                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
2261               
2262                if (p->Valid(mEpsilon))
2263                {
2264                        candidates.push_back(p->CreateReversePolygon());
2265                        DEL_PTR(p);
2266                }
2267        }
2268
2269        // add faces of bounding box (also could be faces of the cell)
2270        for (int i = 0; i < 6; ++ i)
2271        {
2272                VertexContainer vertices;
2273       
2274                for (int j = 0; j < 4; ++ j)
2275                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
2276
2277                candidates.push_back(new Polygon3(vertices));
2278        }
2279
2280        for (int i = 0; i < (int)candidates.size(); ++ i)
2281        {
2282                // polygon is split by all other planes
2283                for (int j = 0; (j < (int)halfSpaces.size()) && candidates[i]; ++ j)
2284                {
2285                        if (i == j) // polygon and plane are coincident
2286                                continue;
2287
2288                        VertexContainer splitPts;
2289                        Polygon3 *frontPoly, *backPoly;
2290
2291                        const int cf = candidates[i]->
2292                                ClassifyPlane(halfSpaces[j], mEpsilon);
2293                       
2294                        switch (cf)
2295                        {
2296                                case Polygon3::SPLIT:
2297                                        frontPoly = new Polygon3();
2298                                        backPoly = new Polygon3();
2299
2300                                        candidates[i]->Split(halfSpaces[j],
2301                                                                                 *frontPoly,
2302                                                                                 *backPoly,
2303                                                                                 mEpsilon);
2304
2305                                        DEL_PTR(candidates[i]);
2306
2307                                        if (frontPoly->Valid(mEpsilon))
2308                                                candidates[i] = frontPoly;
2309                                        else
2310                                                DEL_PTR(frontPoly);
2311
2312                                        DEL_PTR(backPoly);
2313                                        break;
2314                                case Polygon3::BACK_SIDE:
2315                                        DEL_PTR(candidates[i]);
2316                                        break;
2317                                // just take polygon as it is
2318                                case Polygon3::FRONT_SIDE:
2319                                case Polygon3::COINCIDENT:
2320                                default:
2321                                        break;
2322                        }
2323                }
2324               
2325                if (candidates[i])
2326                        geom.mPolys.push_back(candidates[i]);
2327        }
2328}
2329
2330
2331typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;
2332
2333
2334int BspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
2335                                                   const bool onlyUnmailed) const
2336{
2337        stack<bspNodePair> nodeStack;
2338       
2339        BspNodeGeometry nodeGeom;
2340        ConstructGeometry(n, nodeGeom);
2341       
2342        // split planes from the root to this node
2343        // needed to verify that we found neighbor leaf
2344        // TODO: really needed?
2345        vector<Plane3> halfSpaces;
2346        ExtractHalfSpaces(n, halfSpaces);
2347
2348
2349        BspNodeGeometry *rgeom = new BspNodeGeometry();
2350        ConstructGeometry(mRoot, *rgeom);
2351
2352        nodeStack.push(bspNodePair(mRoot, rgeom));
2353
2354        while (!nodeStack.empty())
2355        {
2356                BspNode *node = nodeStack.top().first;
2357                BspNodeGeometry *geom = nodeStack.top().second;
2358       
2359                nodeStack.pop();
2360
2361                if (node->IsLeaf())
2362                {
2363                        // test if this leaf is in valid view space
2364                        if (node->TreeValid() &&
2365                                (node != n) &&
2366                                (!onlyUnmailed || !node->Mailed()))
2367                        {
2368                                bool isAdjacent = true;
2369
2370                                if (1)
2371                                {
2372                                        // test all planes of current node if still adjacent
2373                                        for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
2374                                        {
2375                                                const int cf =
2376                                                        Polygon3::ClassifyPlane(geom->mPolys,
2377                                                                                                        halfSpaces[i],
2378                                                                                                        mEpsilon);
2379
2380                                                if (cf == Polygon3::BACK_SIDE)
2381                                                {
2382                                                        isAdjacent = false;
2383                                                }
2384                                        }
2385                                }
2386                                else
2387                                {
2388                                        // TODO: why is this wrong??
2389                                        // test all planes of current node if still adjacent
2390                                        for (int i = 0; (i < (int)nodeGeom.mPolys.size()) && isAdjacent; ++ i)
2391                                        {
2392                                                Polygon3 *poly = nodeGeom.mPolys[i];
2393
2394                                                const int cf =
2395                                                        Polygon3::ClassifyPlane(geom->mPolys,
2396                                                                                                        poly->GetSupportingPlane(),
2397                                                                                                        mEpsilon);
2398
2399                                                if (cf == Polygon3::BACK_SIDE)
2400                                                {
2401                                                        isAdjacent = false;
2402                                                }
2403                                        }
2404                                }
2405                                // neighbor was found
2406                                if (isAdjacent)
2407                                {       
2408                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2409                                }
2410                        }
2411                }
2412                else
2413                {
2414                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2415
2416                        const int cf = Polygon3::ClassifyPlane(nodeGeom.mPolys,
2417                                                                                                   interior->GetPlane(),
2418                                                                                                   mEpsilon);
2419                       
2420                        BspNode *front = interior->GetFront();
2421                        BspNode *back = interior->GetBack();
2422           
2423                        BspNodeGeometry *fGeom = new BspNodeGeometry();
2424                        BspNodeGeometry *bGeom = new BspNodeGeometry();
2425
2426                        geom->SplitGeometry(*fGeom,
2427                                                                *bGeom,
2428                                                                interior->GetPlane(),
2429                                                                mBox,
2430                                                                mEpsilon);
2431               
2432                        if (cf == Polygon3::FRONT_SIDE)
2433                        {
2434                                nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2435                                DEL_PTR(bGeom);
2436                        }
2437                        else
2438                        {
2439                                if (cf == Polygon3::BACK_SIDE)
2440                                {
2441                                        nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2442                                        DEL_PTR(fGeom);
2443                                }
2444                                else
2445                                {       // random decision
2446                                        nodeStack.push(bspNodePair(front, fGeom));
2447                                        nodeStack.push(bspNodePair(back, bGeom));
2448                                }
2449                        }
2450                }
2451       
2452                DEL_PTR(geom);
2453        }
2454
2455        return (int)neighbors.size();
2456}
2457
2458
2459BspLeaf *BspTree::GetRandomLeaf(const Plane3 &halfspace)
2460{
2461    stack<BspNode *> nodeStack;
2462        nodeStack.push(mRoot);
2463       
2464        int mask = rand();
2465 
2466        while (!nodeStack.empty())
2467        {
2468                BspNode *node = nodeStack.top();
2469                nodeStack.pop();
2470         
2471                if (node->IsLeaf())
2472                {
2473                        return dynamic_cast<BspLeaf *>(node);
2474                }
2475                else
2476                {
2477                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2478                        BspNode *next;
2479       
2480                        BspNodeGeometry geom;
2481                        // todo: not very efficient: constructs full cell everytime
2482                        ConstructGeometry(interior, geom);
2483
2484                        const int cf = Polygon3::ClassifyPlane(geom.mPolys,
2485                                                                                                   halfspace,
2486                                                                                                   mEpsilon);
2487
2488                        if (cf == Polygon3::BACK_SIDE)
2489                                next = interior->GetFront();
2490                        else
2491                                if (cf == Polygon3::FRONT_SIDE)
2492                                        next = interior->GetFront();
2493                        else
2494                        {
2495                                // random decision
2496                                if (mask & 1)
2497                                        next = interior->GetBack();
2498                                else
2499                                        next = interior->GetFront();
2500                                mask = mask >> 1;
2501                        }
2502
2503                        nodeStack.push(next);
2504                }
2505        }
2506       
2507        return NULL;
2508}
2509
2510
2511BspLeaf *BspTree::GetRandomLeaf(const bool onlyUnmailed)
2512{
2513        stack<BspNode *> nodeStack;
2514       
2515        nodeStack.push(mRoot);
2516
2517        int mask = rand();
2518       
2519        while (!nodeStack.empty())
2520        {
2521                BspNode *node = nodeStack.top();
2522                nodeStack.pop();
2523               
2524                if (node->IsLeaf())
2525                {
2526                        if ( (!onlyUnmailed || !node->Mailed()) )
2527                                return dynamic_cast<BspLeaf *>(node);
2528                }
2529                else
2530                {
2531                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
2532
2533                        // random decision
2534                        if (mask & 1)
2535                                nodeStack.push(interior->GetBack());
2536                        else
2537                                nodeStack.push(interior->GetFront());
2538
2539                        mask = mask >> 1;
2540                }
2541        }
2542       
2543        return NULL;
2544}
2545
2546
2547void BspTree::AddToPvs(BspLeaf *leaf,
2548                                           const BoundedRayContainer &rays,
2549                                           int &sampleContributions,
2550                                           int &contributingSamples)
2551{
2552        sampleContributions = 0;
2553        contributingSamples = 0;
2554
2555    BoundedRayContainer::const_iterator it, it_end = rays.end();
2556
2557        ViewCell *vc = leaf->GetViewCell();
2558
2559        // add contributions from samples to the PVS
2560        for (it = rays.begin(); it != it_end; ++ it)
2561        {
2562                int contribution = 0;
2563                Ray *ray = (*it)->mRay;
2564                float relContribution;
2565                if (!ray->intersections.empty())
2566                  contribution += vc->GetPvs().AddSample(ray->intersections[0].mObject,
2567                                                                                                 1.0f,
2568                                                                                                 relContribution);
2569               
2570                if (ray->sourceObject.mObject)
2571                        contribution += vc->GetPvs().AddSample(ray->sourceObject.mObject,
2572                                                                                                   1.0f,
2573                                                                                                   relContribution);
2574               
2575                if (contribution)
2576                {
2577                        sampleContributions += contribution;
2578                        ++ contributingSamples;
2579                }
2580
2581                //if (ray->mFlags & Ray::STORE_BSP_INTERSECTIONS)
2582                //      ray->bspIntersections.push_back(Ray::BspIntersection((*it)->mMinT, this));
2583        }
2584}
2585
2586
2587int BspTree::ComputePvsSize(const BoundedRayContainer &rays) const
2588{
2589        int pvsSize = 0;
2590
2591        BoundedRayContainer::const_iterator rit, rit_end = rays.end();
2592
2593        Intersectable::NewMail();
2594
2595        for (rit = rays.begin(); rit != rays.end(); ++ rit)
2596        {
2597                Ray *ray = (*rit)->mRay;
2598               
2599                if (!ray->intersections.empty())
2600                {
2601                        if (!ray->intersections[0].mObject->Mailed())
2602                        {
2603                                ray->intersections[0].mObject->Mail();
2604                                ++ pvsSize;
2605                        }
2606                }
2607                if (ray->sourceObject.mObject)
2608                {
2609                        if (!ray->sourceObject.mObject->Mailed())
2610                        {
2611                                ray->sourceObject.mObject->Mail();
2612                                ++ pvsSize;
2613                        }
2614                }
2615        }
2616
2617        return pvsSize;
2618}
2619
2620
2621float BspTree::GetEpsilon() const
2622{
2623        return mEpsilon;
2624}
2625
2626
2627int BspTree::CollectMergeCandidates(const vector<BspLeaf *> leaves,
2628                                                                        vector<MergeCandidate> &candidates)
2629{
2630        BspLeaf::NewMail();
2631       
2632        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
2633
2634        int numCandidates = 0;
2635
2636        // find merge candidates and push them into queue
2637        for (it = leaves.begin(); it != it_end; ++ it)
2638        {
2639                BspLeaf *leaf = *it;
2640               
2641                // the same leaves must not be part of two merge candidates
2642                leaf->Mail();
2643                vector<BspLeaf *> neighbors;
2644                FindNeighbors(leaf, neighbors, true);
2645
2646                vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();
2647
2648                // TODO: test if at least one ray goes from one leaf to the other
2649                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
2650                {
2651                        if ((*nit)->GetViewCell() != leaf->GetViewCell())
2652                        {
2653                                MergeCandidate mc(leaf->GetViewCell(), (*nit)->GetViewCell());
2654                                candidates.push_back(mc);
2655
2656                                ++ numCandidates;
2657                                if ((numCandidates % 1000) == 0)
2658                                {
2659                                        cout << "collected " << numCandidates << " merge candidates" << endl;
2660                                }
2661                        }
2662                }
2663        }
2664
2665        Debug << "merge queue: " << (int)candidates.size() << endl;
2666        Debug << "leaves in queue: " << numCandidates << endl;
2667       
2668
2669        return (int)leaves.size();
2670}
2671
2672
2673int BspTree::CollectMergeCandidates(const VssRayContainer &rays,
2674                                                                        vector<MergeCandidate> &candidates)
2675{
2676        ViewCell::NewMail();
2677        long startTime = GetTime();
2678       
2679        map<BspLeaf *, vector<BspLeaf*> > neighborMap;
2680        ViewCellContainer::const_iterator iit;
2681
2682        int numLeaves = 0;
2683       
2684        BspLeaf::NewMail();
2685
2686        for (int i = 0; i < (int)rays.size(); ++ i)
2687        { 
2688                VssRay *ray = rays[i];
2689       
2690                // traverse leaves stored in the rays and compare and
2691                // merge consecutive leaves (i.e., the neighbors in the tree)
2692                if (ray->mViewCells.size() < 2)
2693                        continue;
2694
2695                iit = ray->mViewCells.begin();
2696                BspViewCell *bspVc = dynamic_cast<BspViewCell *>(*(iit ++));
2697                BspLeaf *leaf = bspVc->mLeaf;
2698               
2699                // traverse intersections
2700                // consecutive leaves are neighbors => add them to queue
2701                for (; iit != ray->mViewCells.end(); ++ iit)
2702                {
2703                        // next pair
2704                        BspLeaf *prevLeaf = leaf;
2705                        bspVc = dynamic_cast<BspViewCell *>(*iit);
2706            leaf = bspVc->mLeaf;
2707
2708                        // view space not valid or same view cell
2709                        if (!leaf->TreeValid() || !prevLeaf->TreeValid() ||
2710                                (leaf->GetViewCell() == prevLeaf->GetViewCell()))
2711                                continue;
2712
2713                vector<BspLeaf *> &neighbors = neighborMap[leaf];
2714                       
2715                        bool found = false;
2716
2717                        // both leaves inserted in queue already =>
2718                        // look if double pair already exists
2719                        if (leaf->Mailed() && prevLeaf->Mailed())
2720                        {
2721                                vector<BspLeaf *>::const_iterator it, it_end = neighbors.end();
2722                               
2723                for (it = neighbors.begin(); !found && (it != it_end); ++ it)
2724                                        if (*it == prevLeaf)
2725                                                found = true; // already in queue
2726                        }
2727               
2728                        if (!found)
2729                        {
2730                                // this pair is not in map yet
2731                                // => insert into the neighbor map and the queue
2732                                neighbors.push_back(prevLeaf);
2733                                neighborMap[prevLeaf].push_back(leaf);
2734
2735                                leaf->Mail();
2736                                prevLeaf->Mail();
2737               
2738                                MergeCandidate mc(leaf->GetViewCell(), prevLeaf->GetViewCell());
2739                               
2740                                candidates.push_back(mc);
2741
2742                                if (((int)candidates.size() % 1000) == 0)
2743                                {
2744                                        cout << "collected " << (int)candidates.size() << " merge candidates" << endl;
2745                                }
2746                        }
2747        }
2748        }
2749
2750        Debug << "neighbormap size: " << (int)neighborMap.size() << endl;
2751        Debug << "merge queue: " << (int)candidates.size() << endl;
2752        Debug << "leaves in queue: " << numLeaves << endl;
2753
2754
2755        //-- collect the leaves which haven't been found by ray casting
2756#if 0
2757        cout << "finding additional merge candidates using geometry" << endl;
2758        vector<BspLeaf *> leaves;
2759        CollectLeaves(leaves, true);
2760        Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
2761        CollectMergeCandidates(leaves, candidates);
2762#endif
2763
2764        return numLeaves;
2765}
2766
2767
2768
2769
2770/***************************************************************/
2771/*              BspNodeGeometry Implementation                 */
2772/***************************************************************/
2773
2774
2775BspNodeGeometry::BspNodeGeometry(const BspNodeGeometry &rhs)
2776{
2777        mPolys.reserve(rhs.mPolys.size());
2778       
2779        PolygonContainer::const_iterator it, it_end = rhs.mPolys.end();
2780        for (it = rhs.mPolys.begin(); it != it_end; ++ it)
2781        {
2782                Polygon3 *poly = *it;
2783                mPolys.push_back(new Polygon3(*poly));
2784        }
2785}
2786
2787
2788BspNodeGeometry::~BspNodeGeometry()
2789{
2790        CLEAR_CONTAINER(mPolys);
2791}
2792
2793
2794float BspNodeGeometry::GetArea() const
2795{
2796        return Polygon3::GetArea(mPolys);
2797}
2798
2799
2800float BspNodeGeometry::GetVolume() const
2801{
2802        //-- compute volume using tetrahedralization of the geometry
2803        //   and adding the volume of the single tetrahedrons
2804        float volume = 0;
2805        const float f = 1.0f / 6.0f;
2806
2807        PolygonContainer::const_iterator pit, pit_end = mPolys.end();
2808
2809        // note: can take arbitrary point, e.g., the origin. However,
2810        // we rather take the center of mass to prevents precision errors
2811        const Vector3 center = CenterOfMass();
2812
2813        for (pit = mPolys.begin(); pit != pit_end; ++ pit)
2814        {
2815                Polygon3 *poly = *pit;
2816                const Vector3 v0 = poly->mVertices[0] - center;
2817
2818                for (int i = 1; i < (int)poly->mVertices.size() - 1; ++ i)
2819                {
2820                        const Vector3 v1 = poly->mVertices[i] - center;
2821                        const Vector3 v2 = poly->mVertices[i + 1] - center;
2822
2823                        volume += f * (DotProd(v0, CrossProd(v1, v2)));
2824                }
2825        }
2826
2827        return volume;
2828}
2829
2830
2831Vector3 BspNodeGeometry::CenterOfMass() const
2832{
2833        int n = 0;
2834
2835        Vector3 center(0,0,0);
2836
2837        PolygonContainer::const_iterator pit, pit_end = mPolys.end();
2838
2839        for (pit = mPolys.begin(); pit != pit_end; ++ pit)
2840        {
2841                Polygon3 *poly = *pit;
2842               
2843                VertexContainer::const_iterator vit, vit_end = poly->mVertices.end();
2844
2845                for(vit = poly->mVertices.begin(); vit != vit_end; ++ vit)
2846                {
2847                        center += *vit;
2848                        ++ n;
2849                }
2850        }
2851
2852        return center / (float)n;
2853}
2854
2855
2856void BspNodeGeometry::AddToMesh(Mesh &mesh)
2857{
2858        PolygonContainer::const_iterator it, it_end = mPolys.end();
2859       
2860        for (it = mPolys.begin(); it != mPolys.end(); ++ it)
2861        {
2862                (*it)->AddToMesh(mesh);
2863        }
2864}
2865
2866
2867void BspNodeGeometry::SplitGeometry(BspNodeGeometry &front,
2868                                                                        BspNodeGeometry &back,
2869                                                                        const Plane3 &splitPlane,
2870                                                                        const AxisAlignedBox3 &box,
2871                                                                        const float epsilon) const
2872{       
2873        // get cross section of new polygon
2874        Polygon3 *planePoly = box.CrossSection(splitPlane);
2875
2876        // split polygon with all other polygons
2877        planePoly = SplitPolygon(planePoly, epsilon);
2878
2879        //-- new polygon splits all other polygons
2880        for (int i = 0; i < (int)mPolys.size(); ++ i)
2881        {
2882                /// don't use epsilon here to get exact split planes
2883                const int cf =
2884                        mPolys[i]->ClassifyPlane(splitPlane, Limits::Small);
2885                       
2886                switch (cf)
2887                {
2888                        case Polygon3::SPLIT:
2889                                {
2890                                        Polygon3 *poly = new Polygon3(mPolys[i]->mVertices);
2891
2892                                        Polygon3 *frontPoly = new Polygon3();
2893                                        Polygon3 *backPoly = new Polygon3();
2894                               
2895                                        poly->Split(splitPlane,
2896                                                                *frontPoly,
2897                                                                *backPoly,
2898                                                                epsilon);
2899
2900                                        DEL_PTR(poly);
2901
2902                                        if (frontPoly->Valid(epsilon))
2903                                                front.mPolys.push_back(frontPoly);
2904                                        else
2905                                                DEL_PTR(frontPoly);
2906
2907                                        if (backPoly->Valid(epsilon))
2908                                                back.mPolys.push_back(backPoly);
2909                                        else
2910                                                DEL_PTR(backPoly);
2911                                }
2912                               
2913                                break;
2914                        case Polygon3::BACK_SIDE:
2915                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));                     
2916                                break;
2917                        case Polygon3::FRONT_SIDE:
2918                                front.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));     
2919                                break;
2920                        case Polygon3::COINCIDENT:
2921                                //front.mPolys.push_back(CreateReversePolygon(mPolys[i]));
2922                                back.mPolys.push_back(new Polygon3(mPolys[i]->mVertices));
2923                                break;
2924                        default:
2925                                break;
2926                }
2927        }
2928
2929        //-- finally add the new polygon to the child node geometries
2930        if (planePoly)
2931        {
2932                // add polygon with normal pointing into positive half space to back cell
2933                back.mPolys.push_back(planePoly);
2934                // add polygon with reverse orientation to front cell
2935                front.mPolys.push_back(planePoly->CreateReversePolygon());
2936        }
2937}
2938
2939
2940Polygon3 *BspNodeGeometry::SplitPolygon(Polygon3 *planePoly,
2941                                                                                const float epsilon) const
2942{
2943        if (!planePoly->Valid(epsilon))
2944                DEL_PTR(planePoly);
2945
2946        // polygon is split by all other planes
2947        for (int i = 0; (i < (int)mPolys.size()) && planePoly; ++ i)
2948        {
2949                Plane3 plane = mPolys[i]->GetSupportingPlane();
2950
2951                /// don't use epsilon here to get exact split planes
2952                const int cf =
2953                        planePoly->ClassifyPlane(plane, Limits::Small);
2954                       
2955                // split new polygon with all previous planes
2956                switch (cf)
2957                {
2958                        case Polygon3::SPLIT:
2959                                {
2960                                        Polygon3 *frontPoly = new Polygon3();
2961                                        Polygon3 *backPoly = new Polygon3();
2962
2963                                        planePoly->Split(plane,
2964                                                                         *frontPoly,
2965                                                                         *backPoly,
2966                                                                         epsilon);
2967                                       
2968                                        // don't need anymore
2969                                        DEL_PTR(planePoly);
2970                                        DEL_PTR(frontPoly);
2971
2972                                        // back polygon is belonging to geometry
2973                                        if (backPoly->Valid(epsilon))
2974                                                planePoly = backPoly;
2975                                        else
2976                                                DEL_PTR(backPoly);
2977                                }
2978                                break;
2979                        case Polygon3::FRONT_SIDE:
2980                                DEL_PTR(planePoly);
2981                break;
2982                        // polygon is taken as it is
2983                        case Polygon3::BACK_SIDE:
2984                        case Polygon3::COINCIDENT:
2985                        default:
2986                                break;
2987                }
2988        }
2989
2990        return planePoly;
2991}
2992
2993
2994ViewCell *
2995BspTree::GetViewCell(const Vector3 &point)
2996{
2997  if (mRoot == NULL)
2998        return NULL;
2999 
3000
3001  stack<BspNode *> nodeStack;
3002  nodeStack.push(mRoot);
3003 
3004  ViewCell *viewcell = NULL;
3005 
3006  while (!nodeStack.empty())  {
3007        BspNode *node = nodeStack.top();
3008        nodeStack.pop();
3009       
3010        if (node->IsLeaf()) {
3011          viewcell = dynamic_cast<BspLeaf *>(node)->mViewCell;
3012          break;
3013        } else {
3014         
3015          BspInterior *interior = dynamic_cast<BspInterior *>(node);
3016               
3017          // random decision
3018          if (interior->GetPlane().Side(point) < 0)
3019                nodeStack.push(interior->GetBack());
3020          else
3021                nodeStack.push(interior->GetFront());
3022        }
3023  }
3024 
3025  return viewcell;
3026}
3027
3028
3029void BspNodeGeometry::IncludeInBox(AxisAlignedBox3 &box)
3030{
3031        Polygon3::IncludeInBox(mPolys, box);
3032}
3033
3034
3035bool BspTree::Export(ofstream &stream)
3036{
3037        ExportNode(mRoot, stream);
3038
3039        return true;
3040}
3041
3042
3043void BspTree::ExportNode(BspNode *node, ofstream &stream)
3044{
3045        if (node->IsLeaf())
3046        {
3047                BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
3048                       
3049                int id = -1;
3050                if (leaf->GetViewCell() != mOutOfBoundsCell)
3051                        id = leaf->GetViewCell()->GetId();
3052
3053                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
3054        }
3055        else
3056        {
3057                BspInterior *interior = dynamic_cast<BspInterior *>(node);
3058       
3059                Plane3 plane = interior->GetPlane();
3060                stream << "<Interior plane=\"" << plane.mNormal.x << " "
3061                           << plane.mNormal.y << " " << plane.mNormal.z << " "
3062                           << plane.mD << "\">" << endl;
3063
3064                ExportNode(interior->GetBack(), stream);
3065                ExportNode(interior->GetFront(), stream);
3066
3067                stream << "</Interior>" << endl;
3068        }
3069}
Note: See TracBrowser for help on using the repository browser.