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

Revision 485, 63.5 KB checked in by mattausch, 19 years ago (diff)

changed castlinesegment of vspbpstree
changed rayinfo ray classification for bsp tree
implemented shuffling

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