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

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