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

Revision 589, 69.9 KB checked in by bittner, 18 years ago (diff)

unix2dos on all sources

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