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

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