source: GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellBsp.cpp @ 682

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