source: trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp @ 475

Revision 475, 44.9 KB checked in by mattausch, 19 years ago (diff)

bsp merging possible again

Line 
1#include "Plane3.h"
2#include "VspBspTree.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#include <time.h>
12#include <iomanip>
13#include "Exporter.h"
14#include "Plane3.h"
15#include "ViewCellBsp.h"
16
17//-- static members
18/** Evaluates split plane classification with respect to the plane's
19        contribution for a minimum number of ray splits.
20*/
21const float VspBspTree::sLeastRaySplitsTable[] = {0, 0, 1, 1, 0};
22/** Evaluates split plane classification with respect to the plane's
23        contribution for balanced rays.
24*/
25const float VspBspTree::sBalancedRaysTable[] = {1, -1, 0, 0, 0};
26
27
28int VspBspTree::sFrontId = 0;
29int VspBspTree::sBackId = 0;
30int VspBspTree::sFrontAndBackId = 0;
31
32
33/****************************************************************/
34/*                  class VspBspTree implementation             */
35/****************************************************************/
36
37VspBspTree::VspBspTree():
38mRoot(NULL),
39mPvsUseArea(true),
40mCostNormalizer(Limits::Small)
41{
42        mRootCell = new BspViewCell();
43
44        Randomize(); // initialise random generator for heuristics
45
46        //-- termination criteria for autopartition
47        environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
48        environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
49        environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
50        environment->GetFloatValue("VspBspTree.Termination.minArea", mTermMinArea);     
51        environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
52        environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
53        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
54        environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
55
56        //-- factors for bsp tree split plane heuristics
57        environment->GetFloatValue("VspBspTree.Factor.balancedRays", mBalancedRaysFactor);
58        environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
59        environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
60
61        //-- termination criteria for axis aligned split
62        environment->GetFloatValue("VspBspTree.Termination.AxisAligned.ct_div_ci", mAxisAlignedCtDivCi);
63        environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
64       
65        environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
66                                                         mTermMinRaysForAxisAligned);
67       
68        //-- partition criteria
69        environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
70        environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
71        environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
72        environment->GetFloatValue("VspBspTree.AxisAligned.splitBorder", mAxisAlignedSplitBorder);
73
74        environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
75        environment->GetIntValue("VspBspTree.maxTests", mMaxTests);
76
77        // maximum number of view cells
78        environment->GetIntValue("ViewCells.maxViewCells", mMaxViewCells);
79
80        //--
81        Debug << "******* VSP BSP options ******** " << endl;
82    Debug << "max depth: " << mTermMaxDepth << endl;
83        Debug << "min PVS: " << mTermMinPvs << endl;
84        Debug << "min area: " << mTermMinArea << endl;
85        Debug << "min rays: " << mTermMinRays << endl;
86        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
87        //Debug << "VSP BSP mininam accumulated ray lenght: ", mTermMinAccRayLength) << endl;
88        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
89        Debug << "miss tolerance: " << mTermMissTolerance << endl;
90        Debug << "max view cells: " << mMaxViewCells << endl;
91        Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
92        Debug << "max plane candidates: " << mMaxRayCandidates << endl;
93       
94        Debug << "Split plane strategy: ";
95        if (mSplitPlaneStrategy & RANDOM_POLYGON)
96        {
97                Debug << "random polygon ";
98        }
99        if (mSplitPlaneStrategy & AXIS_ALIGNED)
100        {
101                Debug << "axis aligned ";
102        }
103        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
104        {
105                mCostNormalizer += mLeastRaySplitsFactor;
106                Debug << "least ray splits ";
107        }
108        if (mSplitPlaneStrategy & BALANCED_RAYS)
109        {
110                mCostNormalizer += mBalancedRaysFactor;
111                Debug << "balanced rays ";
112        }
113        if (mSplitPlaneStrategy & PVS)
114        {
115                mCostNormalizer += mPvsFactor;
116                Debug << "pvs";
117        }
118
119        Debug << endl;
120}
121
122
123const BspTreeStatistics &VspBspTree::GetStatistics() const
124{
125        return mStat;
126}
127
128
129VspBspTree::~VspBspTree()
130{
131        DEL_PTR(mRoot);
132        DEL_PTR(mRootCell);
133}
134
135int VspBspTree::AddMeshToPolygons(Mesh *mesh,
136                                                                  PolygonContainer &polys,
137                                                                  MeshInstance *parent)
138{
139        FaceContainer::const_iterator fi;
140       
141        // copy the face data to polygons
142        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
143        {
144                Polygon3 *poly = new Polygon3((*fi), mesh);
145               
146                if (poly->Valid(mEpsilon))
147                {
148                        poly->mParent = parent; // set parent intersectable
149                        polys.push_back(poly);
150                }
151                else
152                        DEL_PTR(poly);
153        }
154        return (int)mesh->mFaces.size();
155}
156
157int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells,
158                                                          PolygonContainer &polys,
159                                                          int maxObjects)
160{
161        int limit = (maxObjects > 0) ?
162                Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
163 
164        int polysSize = 0;
165
166        for (int i = 0; i < limit; ++ i)
167        {
168                if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
169                {
170                        mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
171                        polysSize +=
172                                AddMeshToPolygons(viewCells[i]->GetMesh(),
173                                                                  polys,
174                                                                  viewCells[i]);
175                }
176        }
177
178        return polysSize;
179}
180
181int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects,
182                                                                 PolygonContainer &polys,
183                                                                 int maxObjects)
184{
185        int limit = (maxObjects > 0) ?
186                Min((int)objects.size(), maxObjects) : (int)objects.size();
187 
188        for (int i = 0; i < limit; ++i)
189        {
190                Intersectable *object = objects[i];//*it;
191                Mesh *mesh = NULL;
192
193                switch (object->Type()) // extract the meshes
194                {
195                case Intersectable::MESH_INSTANCE:
196                        mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
197                        break;
198                case Intersectable::VIEW_CELL:
199                        mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
200                        break;
201                        // TODO: handle transformed mesh instances
202                default:
203                        Debug << "intersectable type not supported" << endl;
204                        break;
205                }
206               
207        if (mesh) // copy the mesh data to polygons
208                {
209                        mBox.Include(object->GetBox()); // add to BSP tree aabb
210                        AddMeshToPolygons(mesh, polys, mRootCell);
211                }
212        }
213
214        return (int)polys.size();
215}
216
217void VspBspTree::Construct(const VssRayContainer &sampleRays)
218{
219    mStat.nodes = 1;
220        mBox.Initialize();      // initialise BSP tree bounding box
221       
222        PolygonContainer polys;
223        RayInfoContainer *rays = new RayInfoContainer();
224
225        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
226
227        long startTime = GetTime();
228
229        cout << "Extracting polygons from rays ...";
230
231        Intersectable::NewMail();
232
233        //-- extract polygons intersected by the rays
234        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
235        {
236                VssRay *ray = *rit;
237       
238                if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
239                {
240                        ray->mTerminationObject->Mail();
241                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mTerminationObject);
242                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
243                }
244
245                if (ray->mOriginObject && !ray->mOriginObject->Mailed())
246                {
247                        ray->mOriginObject->Mail();
248                        MeshInstance *obj = dynamic_cast<MeshInstance *>(ray->mOriginObject);
249                        AddMeshToPolygons(obj->GetMesh(), polys, obj);
250                }
251        }
252
253        // compute bounding box
254        Polygon3::IncludeInBox(polys, mBox);
255
256        //-- store rays
257        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
258        {
259                VssRay *ray = *rit;
260               
261                float minT, maxT;
262
263                // TODO: not very efficient to implictly cast between rays types ...
264                if (mBox.GetRaySegment(*ray, minT, maxT))
265                {
266                        float len = ray->Length();
267                       
268                        if (!len)
269                                len = Limits::Small;
270                       
271                        rays->push_back(RayInfo(ray, minT / len, maxT / len));
272                }
273        }
274
275        mTermMinArea *= mBox.SurfaceArea(); // normalize
276        mStat.polys = (int)polys.size();
277
278        cout << "finished" << endl;
279
280        Debug << "\nPolygon extraction: " << (int)polys.size() << " polys extracted from "
281                  << (int)sampleRays.size() << " rays in "
282                  << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;
283
284        Construct(polys, rays);
285
286        // clean up polygons
287        CLEAR_CONTAINER(polys);
288}
289
290void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
291{
292        VspBspTraversalStack tStack;
293
294        mRoot = new BspLeaf();
295
296        // constrruct root node geometry
297        BspNodeGeometry *geom = new BspNodeGeometry();
298        ConstructGeometry(mRoot, *geom);
299
300        VspBspTraversalData tData(mRoot,
301                                                          new PolygonContainer(polys),
302                                                          0,
303                                                          rays,
304                              ComputePvsSize(*rays),
305                                                          geom->GetArea(),
306                                                          geom);
307
308        tStack.push(tData);
309
310        mStat.Start();
311        cout << "Contructing vsp bsp tree ... ";
312
313        long startTime = GetTime();
314       
315        while (!tStack.empty())
316        {
317                tData = tStack.top();
318
319            tStack.pop();
320                // subdivide leaf node
321                BspNode *r = Subdivide(tStack, tData);
322
323                if (r == mRoot)
324                        Debug << "VSP BSP tree construction time spent at root: "
325                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl << endl;
326        }
327
328        cout << "finished\n";
329
330        mStat.Stop();
331}
332
333bool VspBspTree::TerminationCriteriaMet(const VspBspTraversalData &data,
334                                                                                const int numLeaves) const
335{
336        return
337                (((int)data.mRays->size() <= mTermMinRays) ||
338                 (data.mPvs <= mTermMinPvs)   ||
339                 (data.mArea <= mTermMinArea) ||
340                 ((mStat.nodes / 2 + 1) >= mMaxViewCells) ||
341                // (data.GetAvgRayContribution() >= mTermMaxRayContribution) ||
342                 (data.mDepth >= mTermMaxDepth));
343}
344
345BspNode *VspBspTree::Subdivide(VspBspTraversalStack &tStack,
346                                                           VspBspTraversalData &tData)
347{
348        BspNode *newNode = tData.mNode;
349
350        if (!TerminationCriteriaMet(tData, (int)tStack.size()))
351        {
352                PolygonContainer coincident;
353       
354                VspBspTraversalData tFrontData;
355                VspBspTraversalData tBackData;
356
357                // create new interior node and two leaf nodes
358                // or return leaf as it is (if maxCostRatio missed)
359                newNode = SubdivideNode(tData, tFrontData, tBackData, coincident);
360
361                if (!newNode->IsLeaf()) //-- continue subdivision
362                {
363                        // push the children on the stack
364                        tStack.push(tFrontData);
365                        tStack.push(tBackData);
366
367                        // delete old leaf node
368                        DEL_PTR(tData.mNode);   
369                }
370        }
371       
372        //-- terminate traversal 
373        if (newNode->IsLeaf())
374        {
375                BspLeaf *leaf = dynamic_cast<BspLeaf *>(newNode);
376       
377                // create new view cell for this leaf
378                BspViewCell *viewCell = new BspViewCell();
379                leaf->SetViewCell(viewCell);
380                viewCell->mLeaves.push_back(leaf);
381
382                //-- update pvs
383                int conSamp = 0, sampCon = 0;
384                AddToPvs(leaf, *tData.mRays, conSamp, sampCon);
385                       
386                mStat.contributingSamples += conSamp;
387                mStat.sampleContributions += sampCon;
388               
389                EvaluateLeafStats(tData);
390        }
391       
392               
393        //-- cleanup
394        DEL_PTR(tData.mPolygons);
395        DEL_PTR(tData.mRays);
396        DEL_PTR(tData.mGeometry);
397
398        return newNode;
399}
400
401BspNode *VspBspTree::SubdivideNode(VspBspTraversalData &tData,
402                                                                   VspBspTraversalData &frontData,
403                                                                   VspBspTraversalData &backData,
404                                                                   PolygonContainer &coincident)
405{
406        BspLeaf *leaf = dynamic_cast<BspLeaf *>(tData.mNode);
407
408        int maxCostMisses = tData.mMaxCostMisses;
409
410        // select subdivision plane
411        Plane3 splitPlane;
412        if (!SelectPlane(splitPlane, leaf, tData))
413        {
414                ++ maxCostMisses;
415
416                if (maxCostMisses >= mTermMissTolerance)
417                {
418                        // terminate branch because of max cost
419                        ++ mStat.maxCostNodes;
420            return leaf;
421                }
422        }
423
424        mStat.nodes += 2;
425
426        //-- subdivide further
427        BspInterior *interior = new BspInterior(splitPlane);
428
429#ifdef _DEBUG
430        Debug << interior << endl;
431#endif
432       
433        //-- the front and back traversal data is filled with the new values
434        frontData.mPolygons = new PolygonContainer();
435        frontData.mDepth = tData.mDepth + 1;
436        frontData.mRays = new RayInfoContainer();
437        frontData.mGeometry = new BspNodeGeometry();
438
439        backData.mPolygons = new PolygonContainer();
440        backData.mDepth = tData.mDepth + 1;
441        backData.mRays = new RayInfoContainer();
442        backData.mGeometry = new BspNodeGeometry();
443
444        // subdivide rays
445        SplitRays(interior->GetPlane(),
446                          *tData.mRays,
447                          *frontData.mRays,
448                          *backData.mRays);
449       
450        // subdivide polygons
451        mStat.splits += SplitPolygons(interior->GetPlane(),
452                                                                  *tData.mPolygons,
453                                          *frontData.mPolygons,
454                                                                  *backData.mPolygons,
455                                                                  coincident);
456
457       
458        // how often was max cost ratio missed in this branch?
459        frontData.mMaxCostMisses = maxCostMisses;
460        backData.mMaxCostMisses = maxCostMisses;
461
462        // compute pvs
463        frontData.mPvs = ComputePvsSize(*frontData.mRays);
464        backData.mPvs = ComputePvsSize(*backData.mRays);
465
466        // split geometry and compute area
467        if (1)
468        {
469                tData.mGeometry->SplitGeometry(*frontData.mGeometry,
470                                                                           *backData.mGeometry,
471                                                                           interior->GetPlane(),
472                                                                           mBox,
473                                                                           mEpsilon);
474       
475                frontData.mArea = frontData.mGeometry->GetArea();
476                backData.mArea = backData.mGeometry->GetArea();
477        }
478
479        // compute accumulated ray length
480        //frontData.mAccRayLength = AccumulatedRayLength(*frontData.mRays);
481        //backData.mAccRayLength = AccumulatedRayLength(*backData.mRays);
482
483        //-- create front and back leaf
484
485        BspInterior *parent = leaf->GetParent();
486
487        // replace a link from node's parent
488        if (!leaf->IsRoot())
489        {
490                parent->ReplaceChildLink(leaf, interior);
491                interior->SetParent(parent);
492        }
493        else // new root
494        {
495                mRoot = interior;
496        }
497
498        // and setup child links
499        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
500       
501        frontData.mNode = interior->GetFront();
502        backData.mNode = interior->GetBack();
503
504        //DEL_PTR(leaf);
505        return interior;
506}
507
508void VspBspTree::AddToPvs(BspLeaf *leaf,
509                                                  const RayInfoContainer &rays,
510                                                  int &sampleContributions,
511                                                  int &contributingSamples)
512{
513        sampleContributions = 0;
514        contributingSamples = 0;
515
516    RayInfoContainer::const_iterator it, it_end = rays.end();
517
518        ViewCell *vc = leaf->GetViewCell();
519
520        // add contributions from samples to the PVS
521        for (it = rays.begin(); it != it_end; ++ it)
522        {
523                int contribution = 0;
524                VssRay *ray = (*it).mRay;
525                       
526                if (ray->mTerminationObject)
527                        contribution += vc->GetPvs().AddSample(ray->mTerminationObject);
528               
529                if (ray->mOriginObject)
530                        contribution += vc->GetPvs().AddSample(ray->mOriginObject);
531
532                if (contribution)
533                {
534                        sampleContributions += contribution;
535                        ++ contributingSamples;
536                }
537                //leaf->mVssRays.push_back(ray);
538        }
539}
540
541void VspBspTree::SortSplitCandidates(const PolygonContainer &polys,
542                                                                         const int axis,
543                                                                         vector<SortableEntry> &splitCandidates) const
544{
545        splitCandidates.clear();
546 
547        const int requestedSize = 2 * (int)polys.size();
548       
549        // creates a sorted split candidates array 
550        splitCandidates.reserve(requestedSize);
551 
552        PolygonContainer::const_iterator it, it_end = polys.end();
553
554        AxisAlignedBox3 box;
555
556        // insert all queries
557        for(it = polys.begin(); it != it_end; ++ it)
558        {
559                box.Initialize();
560                box.Include(*(*it));
561         
562                splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MIN, box.Min(axis), *it));
563                splitCandidates.push_back(SortableEntry(SortableEntry::POLY_MAX, box.Max(axis), *it));
564        }
565 
566        stable_sort(splitCandidates.begin(), splitCandidates.end());
567}
568
569
570float VspBspTree::BestCostRatio(const PolygonContainer &polys,
571                                                                const AxisAlignedBox3 &box,
572                                                                const int axis,
573                                                                float &position,
574                                                                int &objectsBack,
575                                int &objectsFront) const
576{
577        vector<SortableEntry> splitCandidates;
578
579        SortSplitCandidates(polys, axis, splitCandidates);
580       
581        // go through the lists, count the number of objects left and right
582        // and evaluate the following cost funcion:
583        // C = ct_div_ci  + (ol + or)/queries
584       
585        int objectsLeft = 0, objectsRight = (int)polys.size();
586       
587        float minBox = box.Min(axis);
588        float maxBox = box.Max(axis);
589        float boxArea = box.SurfaceArea();
590 
591        float minBand = minBox + mAxisAlignedSplitBorder * (maxBox - minBox);
592        float maxBand = minBox + (1.0f - mAxisAlignedSplitBorder) * (maxBox - minBox);
593       
594        float minSum = 1e20f;
595        vector<SortableEntry>::const_iterator ci, ci_end = splitCandidates.end();
596
597        for(ci = splitCandidates.begin(); ci != ci_end; ++ ci)
598        {
599                switch ((*ci).type)
600                {
601                        case SortableEntry::POLY_MIN:
602                                ++ objectsLeft;
603                                break;
604                        case SortableEntry::POLY_MAX:
605                            -- objectsRight;
606                                break;
607                        default:
608                                break;
609                }
610               
611                if ((*ci).value > minBand && (*ci).value < maxBand)
612                {
613                        AxisAlignedBox3 lbox = box;
614                        AxisAlignedBox3 rbox = box;
615                        lbox.SetMax(axis, (*ci).value);
616                        rbox.SetMin(axis, (*ci).value);
617
618                        float sum = objectsLeft * lbox.SurfaceArea() +
619                                                objectsRight * rbox.SurfaceArea();
620     
621                        if (sum < minSum)
622                        {
623                                minSum = sum;
624                                position = (*ci).value;
625
626                                objectsBack = objectsLeft;
627                                objectsFront = objectsRight;
628                        }
629                }
630        }
631 
632        const float oldCost = (float)polys.size();
633        const float newCost = mAxisAlignedCtDivCi + minSum / boxArea;
634        const float ratio = newCost / oldCost;
635
636
637#if 0
638  Debug << "====================" << endl;
639  Debug << "costRatio=" << ratio << " pos=" << position<<" t=" << (position - minBox)/(maxBox - minBox)
640      << "\t o=(" << objectsBack << "," << objectsFront << ")" << endl;
641#endif
642  return ratio;
643}
644
645bool VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
646                                        const PolygonContainer &polys) const
647{
648        AxisAlignedBox3 box;
649        box.Initialize();
650       
651        // create bounding box of region
652        Polygon3::IncludeInBox(polys, box);
653       
654        int objectsBack = 0, objectsFront = 0;
655        int axis = 0;
656        float costRatio = MAX_FLOAT;
657        Vector3 position;
658
659        //-- area subdivision
660        for (int i = 0; i < 3; ++ i)
661        {
662                float p = 0;
663                const float r = BestCostRatio(polys, box, i, p, objectsBack, objectsFront);
664               
665                if (r < costRatio)
666                {
667                        costRatio = r;
668                        axis = i;
669                        position = p;
670                }
671        }
672       
673        if (costRatio >= mTermMaxCostRatio)
674                return false;
675
676        Vector3 norm(0,0,0); norm[axis] = 1.0f;
677        plane = Plane3(norm, position);
678
679        return true;
680}
681
682bool VspBspTree::SelectPlane(Plane3 &plane,
683                                                         BspLeaf *leaf,
684                                                         VspBspTraversalData &data)
685{
686        if ((mSplitPlaneStrategy & AXIS_ALIGNED) &&
687                ((int)data.mRays->size() > mTermMinRaysForAxisAligned))
688        {       // TODO: candidates should be rays
689                return SelectAxisAlignedPlane(plane, *data.mPolygons);
690        }
691
692        // simplest strategy: just take next polygon
693        if (mSplitPlaneStrategy & RANDOM_POLYGON)
694        {
695        if (!data.mPolygons->empty())
696                {
697                        const int randIdx = (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
698                        Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
699
700                        plane = nextPoly->GetSupportingPlane();
701                        return true;
702                }
703                else
704                {
705                        //-- choose plane on midpoint of a ray
706                        const int candidateIdx = (int)RandomValue(0, (Real)((int)data.mRays->size() - 1));
707                                                                       
708                        const Vector3 minPt = (*data.mRays)[candidateIdx].ExtrapOrigin();
709                        const Vector3 maxPt = (*data.mRays)[candidateIdx].ExtrapTermination();
710
711                        const Vector3 pt = (maxPt + minPt) * 0.5;
712
713                        const Vector3 normal = (*data.mRays)[candidateIdx].mRay->GetDir();
714                       
715                        plane = Plane3(normal, pt);
716                        return true;
717                }
718
719                return false;
720        }
721
722        // use heuristics to find appropriate plane
723        return SelectPlaneHeuristics(plane, leaf, data);
724}
725
726Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
727{       
728        const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
729       
730        const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
731        const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();
732
733        const Vector3 pt = (maxPt + minPt) * 0.5;
734        const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());
735
736        return Plane3(normal, pt);
737}
738
739Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
740{       
741        Vector3 pt[3];
742       
743        int idx[3];
744        int cmaxT = 0;
745        int cminT = 0;
746        bool chooseMin = false;
747
748        for (int j = 0; j < 3; ++ j)
749        {
750                idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));
751                       
752                if (idx[j] >= (int)rays.size())
753                {
754                        idx[j] -= (int)rays.size();
755                               
756                        chooseMin = (cminT < 2);
757                }
758                else
759                        chooseMin = (cmaxT < 2);
760
761                RayInfo rayInf = rays[idx[j]];
762                pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
763        }       
764
765        return Plane3(pt[0], pt[1], pt[2]);
766}
767
768Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
769{       
770        Vector3 pt[3];
771       
772        int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
773        int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
774
775        // check if rays different
776        if (idx1 == idx2)
777                idx2 = (idx2 + 1) % (int)rays.size();
778
779        const RayInfo ray1 = rays[idx1];
780        const RayInfo ray2 = rays[idx2];
781
782        // normal vector of the plane parallel to both lines
783        const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));
784
785        // vector from line 1 to line 2
786        const Vector3 vd = (ray2.ExtrapOrigin() - ray1.ExtrapOrigin());
787       
788        // project vector on normal to get distance
789        const float dist = DotProd(vd, norm);
790
791        // point on plane lies halfway between the two planes
792        const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;
793
794        return Plane3(norm, planePt);
795}
796
797bool VspBspTree::SelectPlaneHeuristics(Plane3 &bestPlane,
798                                                                           BspLeaf *leaf,
799                                                                           VspBspTraversalData &data)
800{
801        float lowestCost = MAX_FLOAT;
802        // intermediate plane
803        Plane3 plane;
804
805        const int limit = Min((int)data.mPolygons->size(), mMaxPolyCandidates);
806        int maxIdx = (int)data.mPolygons->size();
807       
808        for (int i = 0; i < limit; ++ i)
809        {
810                // assure that no index is taken twice
811                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
812                //Debug << "current Idx: " << maxIdx << " cand idx " << candidateIdx << endl;
813               
814                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
815
816                // swap candidate to the end to avoid testing same plane
817                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
818       
819                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
820
821                // evaluate current candidate
822                const float candidateCost =
823                        SplitPlaneCost(poly->GetSupportingPlane(), data);
824
825                if (candidateCost < lowestCost)
826                {
827                        bestPlane = poly->GetSupportingPlane();
828                        lowestCost = candidateCost;
829                }
830        }
831       
832        //-- choose candidate planes extracted from rays
833        //-- different methods are available
834        for (int i = 0; i < mMaxRayCandidates; ++ i)
835        {
836                plane = ChooseCandidatePlane3(*data.mRays);
837                const float candidateCost = SplitPlaneCost(plane, data);
838
839                if (candidateCost < lowestCost)
840                {
841                        bestPlane = plane;     
842                        lowestCost = candidateCost;
843                }
844        }
845
846#ifdef _DEBUG
847        Debug << "plane lowest cost: " << lowestCost << endl;
848#endif
849
850        // cost ratio miss
851        if (lowestCost > mTermMaxCostRatio)
852                return false;
853
854        return true;
855}
856
857
858inline void VspBspTree::GenerateUniqueIdsForPvs()
859{
860        Intersectable::NewMail(); sBackId = ViewCell::sMailId;
861        Intersectable::NewMail(); sFrontId = ViewCell::sMailId;
862        Intersectable::NewMail(); sFrontAndBackId = ViewCell::sMailId;
863}
864
865float VspBspTree::SplitPlaneCost(const Plane3 &candidatePlane,
866                                                                 const VspBspTraversalData &data)
867{
868        float cost = 0;
869
870        float sumBalancedRays = 0;
871        float sumRaySplits = 0;
872       
873        int frontPvs = 0;
874        int backPvs = 0;
875
876        // probability that view point lies in child
877        float pOverall = 0;
878        float pFront = 0;
879        float pBack = 0;
880
881        const bool pvsUseLen = false;
882
883        if (mSplitPlaneStrategy & PVS)
884        {
885                // create unique ids for pvs heuristics
886                GenerateUniqueIdsForPvs();
887       
888                if (mPvsUseArea) // use front and back cell areas to approximate volume
889                {       
890                        // construct child geometry with regard to the candidate split plane
891                        BspNodeGeometry frontCell;
892                        BspNodeGeometry backCell;
893               
894                        data.mGeometry->SplitGeometry(frontCell,
895                                                                                  backCell,
896                                                                                  candidatePlane,
897                                                                                  mBox,
898                                                                                  mEpsilon);
899               
900                        pFront = frontCell.GetArea();
901                        pBack = backCell.GetArea();
902
903                        pOverall = data.mArea;
904                }
905        }
906               
907        int limit;
908        bool useRand;
909
910        // choose test polyongs randomly if over threshold
911        if ((int)data.mRays->size() > mMaxTests)
912        {
913                useRand = true;
914                limit = mMaxTests;
915        }
916        else
917        {
918                useRand = false;
919                limit = (int)data.mRays->size();
920        }
921
922        for (int i = 0; i < limit; ++ i)
923        {
924                const int testIdx = useRand ? (int)RandomValue(0, (Real)(limit - 1)) : i;
925                RayInfo rayInf = (*data.mRays)[testIdx];
926
927                VssRay *ray = rayInf.mRay;
928                float t;
929                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
930
931                if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
932                {
933                        sumBalancedRays += cf;
934                }
935               
936                if (mSplitPlaneStrategy & BALANCED_RAYS)
937                {
938                        if (cf == 0)
939                                ++ sumRaySplits;
940                }
941
942                if (mSplitPlaneStrategy & PVS)
943                {
944                        // in case the ray intersects an object
945                        // assure that we only count the object
946                        // once for the front and once for the back side of the plane
947                       
948                        // add the termination object
949                        AddObjToPvs(ray->mTerminationObject, cf, frontPvs, backPvs);
950                       
951                        // add the source object
952                        AddObjToPvs(ray->mOriginObject, cf, frontPvs, backPvs);
953                       
954                        // use number or length of rays to approximate volume
955                        if (!mPvsUseArea)
956                        {
957                                float len = 1;
958
959                                if (pvsUseLen) // use length of rays
960                                        len = rayInf.SqrSegmentLength();
961                       
962                                pOverall += len;
963
964                                if (cf == 1)
965                                        pFront += len;
966                                if (cf == -1)
967                                        pBack += len;
968                                if (cf == 0)
969                                {
970                                        // use length of rays to approximate volume
971                                        if (pvsUseLen)
972                                        {
973                                                float newLen = len *
974                                                        (rayInf.GetMaxT() - t) /
975                                                        (rayInf.GetMaxT() - rayInf.GetMinT());
976               
977                                                if (candidatePlane.Side(rayInf.ExtrapOrigin()) <= 0)
978                                                {
979                                                        pBack += newLen;
980                                                        pFront += len - newLen;
981                                                }
982                                                else
983                                                {
984                                                        pFront += newLen;
985                                                        pBack += len - newLen;
986                                                }
987                                        }
988                                        else
989                                        {
990                                                ++ pFront;
991                                                ++ pBack;
992                                        }
993                                }
994                        }
995                }
996        }
997
998        const float raysSize = (float)data.mRays->size() + Limits::Small;
999       
1000        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
1001                cost += mLeastRaySplitsFactor * sumRaySplits / raysSize;
1002
1003        if (mSplitPlaneStrategy & BALANCED_RAYS)
1004                cost += mBalancedRaysFactor * fabs(sumBalancedRays) /  raysSize;
1005       
1006        // pvs criterium
1007        if (mSplitPlaneStrategy & PVS)
1008        {
1009                const float oldCost = pOverall * (float)data.mPvs + Limits::Small;
1010                cost += mPvsFactor * (frontPvs * pFront + backPvs * pBack) / oldCost;
1011
1012                //cost += mPvsFactor * 0.5 * (frontPvs * pFront + backPvs * pBack) / oldCost;
1013                //Debug << "new cost: " << cost << " over" << frontPvs * pFront + backPvs * pBack << " old cost " << oldCost << endl;
1014       
1015                if (0) // give penalty to unbalanced split
1016                        if (((pFront * 0.2 + Limits::Small) > pBack) ||
1017                                (pFront < (pBack * 0.2 + Limits::Small)))
1018                                        cost += 0.5;
1019        }
1020
1021#ifdef _DEBUG
1022        Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
1023                  << " frontpvs: " << frontPvs << " pFront: " << pFront
1024                  << " backpvs: " << backPvs << " pBack: " << pBack << endl << endl;
1025#endif
1026       
1027        // normalize cost by sum of linear factors
1028        return cost / (float)mCostNormalizer;
1029}
1030
1031void VspBspTree::AddObjToPvs(Intersectable *obj,
1032                                                         const int cf,
1033                                                         int &frontPvs,
1034                                                         int &backPvs) const
1035{
1036        if (!obj)
1037                return;
1038        // TODO: does this really belong to no pvs?
1039        //if (cf == Ray::COINCIDENT) return;
1040
1041        // object belongs to both PVS
1042        if (cf >= 0)
1043        {
1044                if ((obj->mMailbox != sFrontId) &&
1045                        (obj->mMailbox != sFrontAndBackId))
1046                {
1047                        ++ frontPvs;
1048
1049                        if (obj->mMailbox == sBackId)
1050                                obj->mMailbox = sFrontAndBackId;       
1051                        else
1052                                obj->mMailbox = sFrontId;                                                               
1053                }
1054        }
1055       
1056        if (cf <= 0)
1057        {
1058                if ((obj->mMailbox != sBackId) &&
1059                        (obj->mMailbox != sFrontAndBackId))
1060                {
1061                        ++ backPvs;
1062
1063                        if (obj->mMailbox == sFrontId)
1064                                obj->mMailbox = sFrontAndBackId;
1065                        else
1066                                obj->mMailbox = sBackId;                               
1067                }
1068        }
1069}
1070
1071void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves) const
1072{
1073        stack<BspNode *> nodeStack;
1074        nodeStack.push(mRoot);
1075 
1076        while (!nodeStack.empty())
1077        {
1078                BspNode *node = nodeStack.top();
1079   
1080                nodeStack.pop();
1081   
1082                if (node->IsLeaf())
1083                {
1084                        BspLeaf *leaf = (BspLeaf *)node;               
1085                        leaves.push_back(leaf);
1086                }
1087                else
1088                {
1089                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1090
1091                        nodeStack.push(interior->GetBack());
1092                        nodeStack.push(interior->GetFront());
1093                }
1094        }
1095}
1096
1097AxisAlignedBox3 VspBspTree::GetBoundingBox() const
1098{
1099        return mBox;
1100}
1101
1102BspNode *VspBspTree::GetRoot() const
1103{
1104        return mRoot;
1105}
1106
1107void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
1108{
1109        // the node became a leaf -> evaluate stats for leafs
1110        BspLeaf *leaf = dynamic_cast<BspLeaf *>(data.mNode);
1111
1112        // store maximal and minimal depth
1113        if (data.mDepth > mStat.maxDepth)
1114                mStat.maxDepth = data.mDepth;
1115
1116        if (data.mDepth < mStat.minDepth)
1117                mStat.minDepth = data.mDepth;
1118       
1119        if (data.mDepth >= mTermMaxDepth)
1120                ++ mStat.maxDepthNodes;
1121
1122        if (data.mPvs < mTermMinPvs)
1123                ++ mStat.minPvsNodes;
1124
1125        if ((int)data.mRays->size() < mTermMinRays)
1126                ++ mStat.minRaysNodes;
1127
1128        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
1129                ++ mStat.maxRayContribNodes;
1130       
1131        if (data.mArea <= mTermMinArea)
1132        {
1133                //Debug << "area: " << data.mArea / mBox.SurfaceArea() << " min area: " << mTermMinArea / mBox.SurfaceArea() << endl;
1134                ++ mStat.minAreaNodes;
1135        }
1136        // accumulate depth to compute average depth
1137        mStat.accumDepth += data.mDepth;
1138
1139#ifdef _DEBUG
1140        Debug << "BSP stats: "
1141                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
1142                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
1143                  << "Area: " << data.mArea << " (min: " << mTermMinArea << "), "
1144                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
1145                  << "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
1146                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
1147#endif
1148}
1149
1150int VspBspTree::CastRay(Ray &ray)
1151{
1152        int hits = 0;
1153 
1154        stack<BspRayTraversalData> tStack;
1155 
1156        float maxt, mint;
1157
1158        if (!mBox.GetRaySegment(ray, mint, maxt))
1159                return 0;
1160
1161        Intersectable::NewMail();
1162
1163        Vector3 entp = ray.Extrap(mint);
1164        Vector3 extp = ray.Extrap(maxt);
1165 
1166        BspNode *node = mRoot;
1167        BspNode *farChild = NULL;
1168       
1169        while (1)
1170        {
1171                if (!node->IsLeaf())
1172                {
1173                        BspInterior *in = dynamic_cast<BspInterior *>(node);
1174
1175                        Plane3 splitPlane = in->GetPlane();
1176                        const int entSide = splitPlane.Side(entp);
1177                        const int extSide = splitPlane.Side(extp);
1178
1179                        if (entSide < 0)
1180                        {
1181                                node = in->GetBack();
1182
1183                                if(extSide <= 0) // plane does not split ray => no far child
1184                                        continue;
1185                                       
1186                                farChild = in->GetFront(); // plane splits ray
1187
1188                        } else if (entSide > 0)
1189                        {
1190                                node = in->GetFront();
1191
1192                                if (extSide >= 0) // plane does not split ray => no far child
1193                                        continue;
1194
1195                                farChild = in->GetBack(); // plane splits ray                   
1196                        }
1197                        else // ray and plane are coincident
1198                        {
1199                                // WHAT TO DO IN THIS CASE ?
1200                                //break;
1201                                node = in->GetFront();
1202                                continue;
1203                        }
1204
1205                        // push data for far child
1206                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
1207
1208                        // find intersection of ray segment with plane
1209                        float t;
1210                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
1211                        maxt *= t;
1212                       
1213                } else // reached leaf => intersection with view cell
1214                {
1215                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1216     
1217                        if (!leaf->GetViewCell()->Mailed())
1218                        {
1219                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
1220                                leaf->GetViewCell()->Mail();
1221                                ++ hits;
1222                        }
1223                       
1224                        //-- fetch the next far child from the stack
1225                        if (tStack.empty())
1226                                break;
1227     
1228                        entp = extp;
1229                        mint = maxt; // NOTE: need this?
1230
1231                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
1232                                break;
1233
1234                        BspRayTraversalData &s = tStack.top();
1235
1236                        node = s.mNode;
1237                        extp = s.mExitPoint;
1238                        maxt = s.mMaxT;
1239
1240                        tStack.pop();
1241                }
1242        }
1243
1244        return hits;
1245}
1246
1247bool VspBspTree::Export(const string filename)
1248{
1249        Exporter *exporter = Exporter::GetExporter(filename);
1250
1251        if (exporter)
1252        {
1253                //exporter->ExportVspBspTree(*this);
1254                return true;
1255        }       
1256
1257        return false;
1258}
1259
1260void VspBspTree::CollectViewCells(ViewCellContainer &viewCells) const
1261{
1262        stack<BspNode *> nodeStack;
1263        nodeStack.push(mRoot);
1264
1265        ViewCell::NewMail();
1266
1267        while (!nodeStack.empty())
1268        {
1269                BspNode *node = nodeStack.top();
1270                nodeStack.pop();
1271
1272                if (node->IsLeaf())
1273                {
1274                        ViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
1275
1276                        if (!viewCell->Mailed())
1277                        {
1278                                viewCell->Mail();
1279                                viewCells.push_back(viewCell);
1280                        }
1281                }
1282                else
1283                {
1284                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1285
1286                        nodeStack.push(interior->GetFront());
1287                        nodeStack.push(interior->GetBack());
1288                }
1289        }
1290}
1291
1292void VspBspTree::EvaluateViewCellsStats(ViewCellsStatistics &stat) const
1293{
1294        stat.Reset();
1295
1296        stack<BspNode *> nodeStack;
1297        nodeStack.push(mRoot);
1298
1299        ViewCell::NewMail();
1300
1301        // exclude root cell
1302        mRootCell->Mail();
1303
1304        while (!nodeStack.empty())
1305        {
1306                BspNode *node = nodeStack.top();
1307                nodeStack.pop();
1308
1309                if (node->IsLeaf())
1310                {
1311                        ++ stat.leaves;
1312
1313                        BspViewCell *viewCell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
1314
1315                        if (!viewCell->Mailed())
1316                        {
1317                                viewCell->Mail();
1318                               
1319                                ++ stat.viewCells;
1320                                const int pvsSize = viewCell->GetPvs().GetSize();
1321
1322                stat.pvs += pvsSize;
1323
1324                                if (pvsSize < 1)
1325                                        ++ stat.emptyPvs;
1326
1327                                if (pvsSize > stat.maxPvs)
1328                                        stat.maxPvs = pvsSize;
1329
1330                                if (pvsSize < stat.minPvs)
1331                                        stat.minPvs = pvsSize;
1332
1333                                if ((int)viewCell->mLeaves.size() > stat.maxLeaves)
1334                                        stat.maxLeaves = (int)viewCell->mLeaves.size();         
1335                        }
1336                }
1337                else
1338                {
1339                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1340
1341                        nodeStack.push(interior->GetFront());
1342                        nodeStack.push(interior->GetBack());
1343                }
1344        }
1345}
1346
1347BspTreeStatistics &VspBspTree::GetStat()
1348{
1349        return mStat;
1350}
1351
1352float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
1353{
1354        float len = 0;
1355
1356        RayInfoContainer::const_iterator it, it_end = rays.end();
1357
1358        for (it = rays.begin(); it != it_end; ++ it)
1359                len += (*it).SegmentLength();
1360
1361        return len;
1362}
1363
1364int VspBspTree::SplitRays(const Plane3 &plane,
1365                                                  RayInfoContainer &rays,
1366                                                  RayInfoContainer &frontRays,
1367                                                  RayInfoContainer &backRays)
1368{
1369        int splits = 0;
1370
1371        while (!rays.empty())
1372        {
1373                RayInfo bRay = rays.back();
1374                rays.pop_back();
1375
1376                VssRay *ray = bRay.mRay;
1377                float t;
1378
1379                        // get classification and receive new t
1380                const int cf = bRay.ComputeRayIntersection(plane, t);
1381       
1382                switch (cf)
1383                {
1384                case -1:
1385                        backRays.push_back(bRay);
1386                        break;
1387                case 1:
1388                        frontRays.push_back(bRay);
1389                        break;
1390                case 0:
1391                        //-- split ray
1392                        //--  look if start point behind or in front of plane
1393                        if (plane.Side(bRay.ExtrapOrigin()) <= 0)
1394                        {       
1395                                backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
1396                                frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
1397                        }
1398                        else
1399                        {
1400                                frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
1401                                backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
1402                        }
1403                        break;
1404                default:
1405                        Debug << "Should not come here 4" << endl;
1406                        break;
1407                }
1408        }
1409
1410        return splits;
1411}
1412
1413void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
1414{
1415        BspNode *lastNode;
1416
1417        do
1418        {
1419                lastNode = n;
1420
1421                // want to get planes defining geometry of this node => don't take
1422                // split plane of node itself
1423                n = n->GetParent();
1424               
1425                if (n)
1426                {
1427                        BspInterior *interior = dynamic_cast<BspInterior *>(n);
1428                        Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
1429
1430            if (interior->GetFront() != lastNode)
1431                                halfSpace.ReverseOrientation();
1432
1433                        halfSpaces.push_back(halfSpace);
1434                }
1435        }
1436        while (n);
1437}
1438
1439void VspBspTree::ConstructGeometry(BspNode *n,
1440                                                                   BspNodeGeometry &cell) const
1441{
1442        PolygonContainer polys;
1443        ConstructGeometry(n, polys);
1444        cell.mPolys = polys;
1445}
1446
1447void VspBspTree::ConstructGeometry(BspNode *n,
1448                                                                   PolygonContainer &cell) const
1449{
1450        vector<Plane3> halfSpaces;
1451        ExtractHalfSpaces(n, halfSpaces);
1452
1453        PolygonContainer candidatePolys;
1454
1455        // bounded planes are added to the polygons (reverse polygons
1456        // as they have to be outfacing
1457        for (int i = 0; i < (int)halfSpaces.size(); ++ i)
1458        {
1459                Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
1460               
1461                if (p->Valid(mEpsilon))
1462                {
1463                        candidatePolys.push_back(p->CreateReversePolygon());
1464                        DEL_PTR(p);
1465                }
1466        }
1467
1468        // add faces of bounding box (also could be faces of the cell)
1469        for (int i = 0; i < 6; ++ i)
1470        {
1471                VertexContainer vertices;
1472       
1473                for (int j = 0; j < 4; ++ j)
1474                        vertices.push_back(mBox.GetFace(i).mVertices[j]);
1475
1476                candidatePolys.push_back(new Polygon3(vertices));
1477        }
1478
1479        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
1480        {
1481                // polygon is split by all other planes
1482                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
1483                {
1484                        if (i == j) // polygon and plane are coincident
1485                                continue;
1486
1487                        VertexContainer splitPts;
1488                        Polygon3 *frontPoly, *backPoly;
1489
1490                        const int cf =
1491                                candidatePolys[i]->ClassifyPlane(halfSpaces[j],
1492                                                                                                 mEpsilon);
1493                       
1494                        switch (cf)
1495                        {
1496                                case Polygon3::SPLIT:
1497                                        frontPoly = new Polygon3();
1498                                        backPoly = new Polygon3();
1499
1500                                        candidatePolys[i]->Split(halfSpaces[j],
1501                                                                                         *frontPoly,
1502                                                                                         *backPoly,
1503                                                                                         mEpsilon);
1504
1505                                        DEL_PTR(candidatePolys[i]);
1506
1507                                        if (frontPoly->Valid(mEpsilon))
1508                                                candidatePolys[i] = frontPoly;
1509                                        else
1510                                                DEL_PTR(frontPoly);
1511
1512                                        DEL_PTR(backPoly);
1513                                        break;
1514                                case Polygon3::BACK_SIDE:
1515                                        DEL_PTR(candidatePolys[i]);
1516                                        break;
1517                                // just take polygon as it is
1518                                case Polygon3::FRONT_SIDE:
1519                                case Polygon3::COINCIDENT:
1520                                default:
1521                                        break;
1522                        }
1523                }
1524               
1525                if (candidatePolys[i])
1526                        cell.push_back(candidatePolys[i]);
1527        }
1528}
1529
1530void VspBspTree::ConstructGeometry(BspViewCell *vc, PolygonContainer &vcGeom) const
1531{
1532        vector<BspLeaf *> leaves = vc->mLeaves;
1533
1534        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
1535
1536        for (it = leaves.begin(); it != it_end; ++ it)
1537                ConstructGeometry(*it, vcGeom);
1538}
1539
1540int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
1541                                                   const bool onlyUnmailed) const
1542{
1543        PolygonContainer cell;
1544
1545        ConstructGeometry(n, cell);
1546
1547        stack<BspNode *> nodeStack;
1548        nodeStack.push(mRoot);
1549               
1550        // planes needed to verify that we found neighbor leaf.
1551        vector<Plane3> halfSpaces;
1552        ExtractHalfSpaces(n, halfSpaces);
1553
1554        while (!nodeStack.empty())
1555        {
1556                BspNode *node = nodeStack.top();
1557                nodeStack.pop();
1558
1559                if (node->IsLeaf())
1560                {
1561            if (node != n && (!onlyUnmailed || !node->Mailed()))
1562                        {
1563                                // test all planes of current node if candidate really
1564                                // is neighbour
1565                                PolygonContainer neighborCandidate;
1566                                ConstructGeometry(node, neighborCandidate);
1567                               
1568                                bool isAdjacent = true;
1569                                for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
1570                                {
1571                                        const int cf =
1572                                                Polygon3::ClassifyPlane(neighborCandidate,
1573                                                                                                halfSpaces[i],
1574                                                                                                mEpsilon);
1575
1576                                        if (cf == Polygon3::BACK_SIDE)
1577                                                isAdjacent = false;
1578                                }
1579
1580                                if (isAdjacent)
1581                                        neighbors.push_back(dynamic_cast<BspLeaf *>(node));
1582
1583                                CLEAR_CONTAINER(neighborCandidate);
1584                        }
1585                }
1586                else
1587                {
1588                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1589       
1590                        const int cf = Polygon3::ClassifyPlane(cell,
1591                                                                                                   interior->GetPlane(),
1592                                                                                                   mEpsilon);
1593
1594                        if (cf == Polygon3::FRONT_SIDE)
1595                                nodeStack.push(interior->GetFront());
1596                        else
1597                                if (cf == Polygon3::BACK_SIDE)
1598                                        nodeStack.push(interior->GetBack());
1599                                else
1600                                {
1601                                        // random decision
1602                                        nodeStack.push(interior->GetBack());
1603                                        nodeStack.push(interior->GetFront());
1604                                }
1605                }
1606        }
1607       
1608        CLEAR_CONTAINER(cell);
1609        return (int)neighbors.size();
1610}
1611
1612BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
1613{
1614    stack<BspNode *> nodeStack;
1615        nodeStack.push(mRoot);
1616       
1617        int mask = rand();
1618 
1619        while (!nodeStack.empty())
1620        {
1621                BspNode *node = nodeStack.top();
1622                nodeStack.pop();
1623         
1624                if (node->IsLeaf())
1625                {
1626                        return dynamic_cast<BspLeaf *>(node);
1627                }
1628                else
1629                {
1630                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1631                       
1632                        BspNode *next;
1633       
1634                        PolygonContainer cell;
1635
1636                        // todo: not very efficient: constructs full cell everytime
1637                        ConstructGeometry(interior, cell);
1638
1639                        const int cf = Polygon3::ClassifyPlane(cell, halfspace, mEpsilon);
1640
1641                        if (cf == Polygon3::BACK_SIDE)
1642                                next = interior->GetFront();
1643                        else
1644                                if (cf == Polygon3::FRONT_SIDE)
1645                                        next = interior->GetFront();
1646                        else
1647                        {
1648                                // random decision
1649                                if (mask & 1)
1650                                        next = interior->GetBack();
1651                                else
1652                                        next = interior->GetFront();
1653                                mask = mask >> 1;
1654                        }
1655
1656                        nodeStack.push(next);
1657                }
1658        }
1659       
1660        return NULL;
1661}
1662
1663BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
1664{
1665        stack<BspNode *> nodeStack;
1666       
1667        nodeStack.push(mRoot);
1668
1669        int mask = rand();
1670       
1671        while (!nodeStack.empty())
1672        {
1673                BspNode *node = nodeStack.top();
1674                nodeStack.pop();
1675               
1676                if (node->IsLeaf())
1677                {
1678                        if ( (!onlyUnmailed || !node->Mailed()) )
1679                                return dynamic_cast<BspLeaf *>(node);
1680                }
1681                else
1682                {
1683                        BspInterior *interior = dynamic_cast<BspInterior *>(node);
1684
1685                        // random decision
1686                        if (mask & 1)
1687                                nodeStack.push(interior->GetBack());
1688                        else
1689                                nodeStack.push(interior->GetFront());
1690
1691                        mask = mask >> 1;
1692                }
1693        }
1694       
1695        return NULL;
1696}
1697
1698int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
1699{
1700        int pvsSize = 0;
1701
1702        RayInfoContainer::const_iterator rit, rit_end = rays.end();
1703
1704        Intersectable::NewMail();
1705
1706        for (rit = rays.begin(); rit != rays.end(); ++ rit)
1707        {
1708                VssRay *ray = (*rit).mRay;
1709               
1710                if (ray->mOriginObject)
1711                {
1712                        if (!ray->mOriginObject->Mailed())
1713                        {
1714                                ray->mOriginObject->Mail();
1715                                ++ pvsSize;
1716                        }
1717                }
1718                if (ray->mTerminationObject)
1719                {
1720                        if (!ray->mTerminationObject->Mailed())
1721                        {
1722                                ray->mTerminationObject->Mail();
1723                                ++ pvsSize;
1724                        }
1725                }
1726        }
1727
1728        return pvsSize;
1729}
1730
1731float VspBspTree::GetEpsilon() const
1732{
1733        return mEpsilon;
1734}
1735
1736BspViewCell *VspBspTree::GetRootCell() const
1737{
1738        return mRootCell;
1739}
1740
1741int VspBspTree::SplitPolygons(const Plane3 &plane,
1742                                                          PolygonContainer &polys,
1743                                                          PolygonContainer &frontPolys,
1744                                                          PolygonContainer &backPolys,
1745                                                          PolygonContainer &coincident) const
1746{
1747        int splits = 0;
1748
1749        while (!polys.empty())
1750        {
1751                Polygon3 *poly = polys.back();
1752                polys.pop_back();
1753
1754                // classify polygon
1755                const int cf = poly->ClassifyPlane(plane, mEpsilon);
1756
1757                switch (cf)
1758                {
1759                        case Polygon3::COINCIDENT:
1760                                coincident.push_back(poly);
1761                                break;                 
1762                        case Polygon3::FRONT_SIDE:     
1763                                frontPolys.push_back(poly);
1764                                break;
1765                        case Polygon3::BACK_SIDE:
1766                                backPolys.push_back(poly);
1767                                break;
1768                        case Polygon3::SPLIT:
1769                                backPolys.push_back(poly);
1770                                frontPolys.push_back(poly);
1771                                ++ splits;
1772                                break;
1773                        default:
1774                Debug << "SHOULD NEVER COME HERE\n";
1775                                break;
1776                }
1777        }
1778
1779        return splits;
1780}
1781
1782
1783int VspBspTree::CastLineSegment(const Vector3 &origin,
1784                                                                const Vector3 &termination,
1785                                                                vector<ViewCell *> &viewcells)
1786{
1787        int hits = 0;
1788        stack<BspRayTraversalData> tStack;
1789 
1790        float mint = 0.0f, maxt = 1.0f;
1791 
1792        Intersectable::NewMail();
1793 
1794        Vector3 entp = origin;
1795        Vector3 extp = termination;
1796 
1797        BspNode *node = mRoot;
1798        BspNode *farChild = NULL;
1799 
1800        while (1)
1801        {
1802                if (!node->IsLeaf()) 
1803                {
1804                        BspInterior *in = dynamic_cast<BspInterior *>(node);
1805         
1806                        Plane3 splitPlane = in->GetPlane();
1807                        const int entSide = splitPlane.Side(entp);
1808                        const int extSide = splitPlane.Side(extp);
1809         
1810                        if (entSide < 0)
1811                        {
1812                                node = in->GetBack();
1813               
1814                                if(extSide <= 0) // plane does not split ray => no far child
1815                                        continue;
1816               
1817                                farChild = in->GetFront(); // plane splits ray
1818                        } else if (entSide > 0)
1819                        {
1820                                node = in->GetFront();
1821                 
1822                                if (extSide >= 0) // plane does not split ray => no far child
1823                                        continue;
1824                 
1825                                farChild = in->GetBack(); // plane splits ray                   
1826                        }
1827                        else // ray and plane are coincident
1828                        {
1829                                // WHAT TO DO IN THIS CASE ?
1830                                //break;
1831                                node = in->GetFront();
1832                                continue;
1833                        }
1834         
1835                        // push data for far child
1836                        tStack.push(BspRayTraversalData(farChild, extp, maxt));
1837         
1838                        // find intersection of ray segment with plane
1839                        float t;
1840                        extp = splitPlane.FindIntersection(origin, extp, &t);
1841                        maxt *= t;
1842         
1843                } else
1844                {
1845                        // reached leaf => intersection with view cell
1846                        BspLeaf *leaf = dynamic_cast<BspLeaf *>(node);
1847         
1848                        if (!leaf->GetViewCell()->Mailed())
1849                        {
1850                                viewcells.push_back(leaf->GetViewCell());
1851                                leaf->GetViewCell()->Mail();
1852                                ++ hits;
1853                        }
1854         
1855                        //-- fetch the next far child from the stack
1856                        if (tStack.empty())
1857                                break;
1858     
1859                        entp = extp;
1860                        mint = maxt; // NOTE: need this?
1861         
1862                       
1863                        BspRayTraversalData &s = tStack.top();
1864           
1865                        node = s.mNode;
1866                        extp = s.mExitPoint;
1867                        maxt = s.mMaxT;
1868         
1869                        tStack.pop();
1870                }
1871        }
1872        return hits;
1873}
Note: See TracBrowser for help on using the repository browser.