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

Revision 495, 67.9 KB checked in by mattausch, 19 years ago (diff)

fixed bug in tree collapse

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