source: GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.cpp @ 652

Revision 652, 77.4 KB checked in by mattausch, 19 years ago (diff)

started to implement split plane queue

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