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

Revision 2071, 98.6 KB checked in by mattausch, 17 years ago (diff)
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"
[1603]19#include "IntersectableWrapper.h"
[463]20
[1077]21
22
[863]23namespace GtpVisibilityPreprocessor {
[860]24
[1077]25
[639]26#define USE_FIXEDPOINT_T 0
[1020]27#define COUNT_ORIGIN_OBJECTS 1
[590]28
[2071]29#define STORE_PVS 1
[1580]30
[2066]31
[1580]32//////////////
[463]33//-- static members
[508]34
[482]35int VspBspTree::sFrontId = 0;
[463]36int VspBspTree::sBackId = 0;
37int VspBspTree::sFrontAndBackId = 0;
38
[1020]39
40
[697]41typedef pair<BspNode *, BspNodeGeometry *> bspNodePair;
[463]42
[557]43
[579]44// pvs penalty can be different from pvs size
[1765]45inline static float EvalPvsPenalty(const float pvs,
46                                                                   const float lower,
47                                                                   const float upper)
[579]48{
49        // clamp to minmax values
50        if (pvs < lower)
51                return (float)lower;
52        if (pvs > upper)
53                return (float)upper;
[578]54
[579]55        return (float)pvs;
56}
[463]57
[557]58
[579]59
60
[601]61/******************************************************************************/
62/*                       class VspBspTree implementation                      */
63/******************************************************************************/
[579]64
65
[1004]66VspBspTree::VspBspTree():
[463]67mRoot(NULL),
[547]68mUseAreaForPvs(false),
[478]69mCostNormalizer(Limits::Small),
70mViewCellsManager(NULL),
[497]71mOutOfBoundsCell(NULL),
[643]72mStoreRays(false),
[605]73mRenderCostWeight(0.5),
[610]74mUseRandomAxis(false),
[650]75mTimeStamp(1)
[463]76{
[486]77        bool randomize = false;
[1004]78        Environment::GetSingleton()->GetBoolValue("VspBspTree.Construction.randomize", randomize);
[1011]79        if (randomize) Randomize(); // initialise random generator for heuristics
[463]80
[1557]81        //////////////////
[463]82        //-- termination criteria for autopartition
[1557]83
[1004]84        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
85        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
86        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
87        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.minProbability", mTermMinProbability);
88        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
89        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
[1023]90       
[1004]91        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
92        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
[666]93
[1557]94        ////////////////////////
95        //-- cost ratios for early tree termination
[1004]96        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
97        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
98        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
[610]99
[1557]100        ///////////
[463]101        //-- factors for bsp tree split plane heuristics
[1557]102
[1004]103        Environment::GetSingleton()->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
104        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
[463]105
[1557]106        //////////
107        //-- partition criteria
[482]108
[1004]109        Environment::GetSingleton()->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
110        Environment::GetSingleton()->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
111        Environment::GetSingleton()->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
[463]112
[1004]113        Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
114        Environment::GetSingleton()->GetIntValue("VspBspTree.maxTests", mMaxTests);
[463]115
[564]116        // if only the driving axis is used for axis aligned split
[1004]117        Environment::GetSingleton()->GetBoolValue("VspBspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
[580]118       
[1557]119        //////////////////////
[508]120        //-- termination criteria for axis aligned split
[1004]121        Environment::GetSingleton()->GetFloatValue("VspBspTree.Termination.AxisAligned.maxRayContribution",
[665]122                                                                mTermMaxRayContriForAxisAligned);
[1004]123        Environment::GetSingleton()->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
[665]124                                                         mTermMinRaysForAxisAligned);
[487]125
[1004]126        Environment::GetSingleton()->GetFloatValue("VspBspTree.maxStaticMemory", mMaxMemory);
[508]127
[1004]128        Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.renderCostWeight", mRenderCostWeight);
[1020]129        Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.renderCostDecreaseWeight", mRenderCostDecreaseWeight);
[1004]130        Environment::GetSingleton()->GetBoolValue("VspBspTree.usePolygonSplitIfAvailable", mUsePolygonSplitIfAvailable);
[547]131
[1004]132        Environment::GetSingleton()->GetBoolValue("VspBspTree.useCostHeuristics", mUseCostHeuristics);
133        Environment::GetSingleton()->GetBoolValue("VspBspTree.useSplitCostQueue", mUseSplitCostQueue);
134        Environment::GetSingleton()->GetBoolValue("VspBspTree.simulateOctree", mCirculatingAxis);
135        Environment::GetSingleton()->GetBoolValue("VspBspTree.useRandomAxis", mUseRandomAxis);
136        Environment::GetSingleton()->GetIntValue("VspBspTree.nodePriorityQueueType", mNodePriorityQueueType);
[580]137
[676]138       
[660]139        char subdivisionStatsLog[100];
[1004]140        Environment::GetSingleton()->GetStringValue("VspBspTree.subdivisionStats", subdivisionStatsLog);
[660]141        mSubdivisionStats.open(subdivisionStatsLog);
[612]142
[1004]143        Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.minBand", mMinBand);
144        Environment::GetSingleton()->GetFloatValue("VspBspTree.Construction.maxBand", mMaxBand);
[1006]145        Environment::GetSingleton()->GetBoolValue("VspBspTree.Construction.useDrivingAxisForMaxCost", mUseDrivingAxisIfMaxCostViolated);
[801]146
[1557]147        /////////
[478]148        //-- debug output
[580]149
[473]150        Debug << "******* VSP BSP options ******** " << endl;
151    Debug << "max depth: " << mTermMaxDepth << endl;
152        Debug << "min PVS: " << mTermMinPvs << endl;
[547]153        Debug << "min probabiliy: " << mTermMinProbability << endl;
[473]154        Debug << "min rays: " << mTermMinRays << endl;
155        Debug << "max ray contri: " << mTermMaxRayContribution << endl;
156        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
157        Debug << "miss tolerance: " << mTermMissTolerance << endl;
158        Debug << "max view cells: " << mMaxViewCells << endl;
159        Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
[486]160        Debug << "randomize: " << randomize << endl;
[582]161
[551]162        Debug << "using area for pvs: " << mUseAreaForPvs << endl;
[580]163        Debug << "render cost weight: " << mRenderCostWeight << endl;
[663]164        Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
165        Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
166        Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
167        Debug << "max memory: " << mMaxMemory << endl;
168        Debug << "use poly split if available: " << mUsePolygonSplitIfAvailable << endl;
169        Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
170        Debug << "use split cost queue: " << mUseSplitCostQueue << endl;
171        Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
[664]172        Debug << "use random axis: " << mUseRandomAxis << endl;
[735]173        Debug << "priority queue type: " << mNodePriorityQueueType << endl;
[822]174        Debug << "circulating axis: " << mCirculatingAxis << endl;
[801]175        Debug << "minband: " << mMinBand << endl;
176        Debug << "maxband: " << mMaxBand << endl;
[1006]177        Debug << "use driving axis for max cost: " << mUseDrivingAxisIfMaxCostViolated << endl;
[1020]178        Debug << "render cost decrease weight: " << mRenderCostDecreaseWeight << endl;
[663]179
[463]180        Debug << "Split plane strategy: ";
[564]181
[463]182        if (mSplitPlaneStrategy & RANDOM_POLYGON)
[474]183        {
[463]184                Debug << "random polygon ";
[474]185        }
[463]186        if (mSplitPlaneStrategy & AXIS_ALIGNED)
[472]187        {
[463]188                Debug << "axis aligned ";
[472]189        }
[665]190        if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
[472]191        {
[474]192                mCostNormalizer += mLeastRaySplitsFactor;
[463]193                Debug << "least ray splits ";
[472]194        }
[463]195        if (mSplitPlaneStrategy & BALANCED_RAYS)
[472]196        {
[474]197                mCostNormalizer += mBalancedRaysFactor;
[463]198                Debug << "balanced rays ";
[472]199        }
[463]200        if (mSplitPlaneStrategy & PVS)
[472]201        {
[474]202                mCostNormalizer += mPvsFactor;
[463]203                Debug << "pvs";
[665]204        }
[482]205
[1557]206        Debug << endl;
[489]207
[1233]208        mLocalSubdivisionCandidates = new vector<SortableEntry>;
[463]209}
[580]210
[1027]211
[508]212BspViewCell *VspBspTree::GetOutOfBoundsCell()
213{
214        return mOutOfBoundsCell;
215}
[463]216
[508]217
[489]218BspViewCell *VspBspTree::GetOrCreateOutOfBoundsCell()
219{
220        if (!mOutOfBoundsCell)
[508]221        {
[489]222                mOutOfBoundsCell = new BspViewCell();
[1557]223                mOutOfBoundsCell->SetId(OUT_OF_BOUNDS_ID);
[547]224                mOutOfBoundsCell->SetValid(false);
[508]225        }
[547]226
[489]227        return mOutOfBoundsCell;
228}
229
230
[482]231const BspTreeStatistics &VspBspTree::GetStatistics() const
[463]232{
[574]233        return mBspStats;
[463]234}
235
236
237VspBspTree::~VspBspTree()
238{
239        DEL_PTR(mRoot);
[1233]240        DEL_PTR(mLocalSubdivisionCandidates);
[463]241}
242
[579]243
[482]244int VspBspTree::AddMeshToPolygons(Mesh *mesh,
[1632]245                                                                  PolygonContainer &polys) const
[463]246{
[1632]247        if (!mesh) return 0;
[1603]248
[463]249        FaceContainer::const_iterator fi;
[482]250
[463]251        // copy the face data to polygons
252        for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
253        {
254                Polygon3 *poly = new Polygon3((*fi), mesh);
[482]255
[463]256                if (poly->Valid(mEpsilon))
257                {
258                        polys.push_back(poly);
259                }
260                else
[1558]261                {
[463]262                        DEL_PTR(poly);
[1558]263                }
[463]264        }
[1603]265
[463]266        return (int)mesh->mFaces.size();
267}
268
[580]269
[1632]270void VspBspTree::ExtractPolygons(Intersectable *object, PolygonContainer &polys) const
[463]271{
[1632]272        // extract the polygons from the intersectables
273        switch (object->Type())
[463]274        {
275                case Intersectable::MESH_INSTANCE:
[1632]276                        {
[2017]277                                Mesh *mesh = static_cast<MeshInstance *>(object)->GetMesh();
[1632]278                                AddMeshToPolygons(mesh, polys);
279                        }
[463]280                        break;
281                case Intersectable::VIEW_CELL:
[1632]282                        {
[2017]283                                Mesh *mesh = static_cast<ViewCell *>(object)->GetMesh();
[1632]284                                AddMeshToPolygons(mesh, polys);
285                                break;
286                        }
[1001]287                case Intersectable::TRANSFORMED_MESH_INSTANCE:
288                        {
[1603]289                                TransformedMeshInstance *mi =
[2017]290                                        static_cast<TransformedMeshInstance *>(object);
[1001]291
[1603]292                                Mesh mesh;
293                                mi->GetTransformedMesh(mesh);
[1632]294                                AddMeshToPolygons(&mesh, polys);
[1001]295                        }
[1632]296                        break;
[1603]297                case Intersectable::TRIANGLE_INTERSECTABLE:
298                        {
299                                TriangleIntersectable *intersect =
[2017]300                                        static_cast<TriangleIntersectable *>(object);
[1603]301
302                                Polygon3 *poly = new Polygon3(intersect->GetItem());
303
304                                if (poly->Valid(mEpsilon))     
[1632]305                                {
[1603]306                                        polys.push_back(poly);
[1632]307                                }
[1603]308                                else
[1632]309                                {
[1603]310                                        delete poly;
[1632]311                                }
[1603]312                        }
[1632]313                        break;
[463]314                default:
[1632]315                        Debug << "intersectable type not supported" << endl;
[463]316                        break;
[1632]317        }
318}
[482]319
[1632]320
321int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects,
322                                                                 PolygonContainer &polys,
323                                                                 int maxObjects)
324{
325        const int limit = (maxObjects > 0) ?
326                Min((int)objects.size(), maxObjects) : (int)objects.size();
327
328        for (int i = 0; i < limit; ++i)
329        {
330                Intersectable *object = objects[i];//*it;
331                ExtractPolygons(object, polys);
332
[2070]333                 // add to BSP tree aabb
334                mBoundingBox.Include(object->GetBox());
[463]335        }
336
337        return (int)polys.size();
338}
339
[580]340
[1027]341void VspBspTree::ComputeBoundingBox(const VssRayContainer &sampleRays,
342                                                                        AxisAlignedBox3 *forcedBoundingBox)
343{
344        if (forcedBoundingBox)
345        {
[1563]346                mBoundingBox = *forcedBoundingBox;
[1027]347        }
348        else // compute vsp tree bounding box
349        {
[1563]350                mBoundingBox.Initialize();
[1027]351
352                VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
353
354                //-- compute bounding box
355        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
356                {
357                        VssRay *ray = *rit;
358
359                        // compute bounding box of view space
[1563]360                        mBoundingBox.Include(ray->GetTermination());
361                        mBoundingBox.Include(ray->GetOrigin());
[1027]362                }
363        }
364}
365
366
[483]367void VspBspTree::Construct(const VssRayContainer &sampleRays,
368                                                   AxisAlignedBox3 *forcedBoundingBox)
[463]369{
[1027]370        // Compute the bounding box from the rays
371        ComputeBoundingBox(sampleRays, forcedBoundingBox);
[484]372       
[463]373        PolygonContainer polys;
374        RayInfoContainer *rays = new RayInfoContainer();
375
[1564]376        ////////////
377        //-- extract polygons from rays if polygon candidate planes are required
378
[1027]379        if (mMaxPolyCandidates)
380        {
381                int numObj = 0;
382                Intersectable::NewMail();
[463]383
[1027]384        cout << "Extracting polygons from rays ... ";
[1580]385                const long startTime = GetTime();
[463]386
[1027]387        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
[542]388
[1027]389                //-- extract polygons intersected by the rays
390                for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
[463]391                {
[1027]392                        VssRay *ray = *rit;
393                        Intersectable *obj = ray->mTerminationObject;
[1002]394
[1632]395                        ++ numObj;
[1027]396
[1632]397                        /////////
398                        //-- compute bounding box
[1603]399
[1632]400                        if (!forcedBoundingBox)
401                        {
402                                mBoundingBox.Include(ray->mTermination);
[1002]403                        }
[1027]404
[1563]405                        if ((mBoundingBox.IsInside(ray->mOrigin) || !forcedBoundingBox) &&
[1027]406                                ray->mOriginObject &&
407                                !ray->mOriginObject->Mailed())
408                        {               
409                                ray->mOriginObject->Mail();
[1632]410                                ExtractPolygons(ray->mOriginObject, polys);
411                                                               
[1027]412                                ++ numObj;
[1002]413                        }
[1027]414                }
[1002]415
[1027]416                // throw out unnecessary polygons
417                PreprocessPolygons(polys);
418                cout << "finished" << endl;
[463]419
[1027]420                Debug << "\n" << (int)polys.size() << " polys extracted from "
421                  << (int)sampleRays.size() << " rays in "
422                  << TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;
[463]423        }
[535]424       
[1027]425        Debug << "maximal pvs (i.e., pvs still considered as valid): "
[653]426                  << mViewCellsManager->GetMaxPvsSize() << endl;
[580]427
[1027]428        VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
429
[2070]430        /////////
[463]431        //-- store rays
432        for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
433        {
434                VssRay *ray = *rit;
[482]435
[463]436                float minT, maxT;
437
[564]438                static Ray hray;
439                hray.Init(*ray);
[562]440
[483]441                // TODO: not very efficient to implictly cast between rays types
[1563]442                if (mBoundingBox.GetRaySegment(hray, minT, maxT))
[463]443                {
444                        float len = ray->Length();
[482]445
446                        if (!len)
[463]447                                len = Limits::Small;
[482]448
[463]449                        rays->push_back(RayInfo(ray, minT / len, maxT / len));
450                }
451        }
452
[587]453        // normalize
[547]454        if (mUseAreaForPvs)
[1563]455                mTermMinProbability *= mBoundingBox.SurfaceArea();
[547]456        else
[1563]457                mTermMinProbability *= mBoundingBox.GetVolume();
[547]458
[648]459
[1027]460        mBspStats.nodes = 1;
[574]461        mBspStats.polys = (int)polys.size();
[1449]462        mBspStats.mGlobalCostMisses = 0;
[654]463
[463]464
[655]465        // use split cost priority queue
466        if (mUseSplitCostQueue)
467        {
[654]468                ConstructWithSplitQueue(polys, rays);
[655]469        }
[654]470        else
[655]471        {
[654]472                Construct(polys, rays);
[655]473        }
[463]474
475        // clean up polygons
476        CLEAR_CONTAINER(polys);
477}
478
[508]479
[612]480// TODO: return memory usage in MB
[656]481float VspBspTree::GetMemUsage() const
[508]482{
[656]483        return (float)
484                 (sizeof(VspBspTree) +
485                  mBspStats.Leaves() * sizeof(BspLeaf) +
486                  mCreatedViewCells * sizeof(BspViewCell) +
[1189]487                  mBspStats.pvs * sizeof(PvsData) +
[656]488                  mBspStats.Interior() * sizeof(BspInterior) +
489                  mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
[508]490}
491
492
[463]493void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
494{
[600]495        VspBspTraversalQueue tQueue;
[463]496
[1016]497        /// create new vsp tree
[463]498        mRoot = new BspLeaf();
499
500        // constrruct root node geometry
501        BspNodeGeometry *geom = new BspNodeGeometry();
502        ConstructGeometry(mRoot, *geom);
503
[1016]504        /// we use the overall probability as normalizer
505        /// either the overall area or the volume
[547]506        const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();
507
[1020]508        /// first traversal data
[482]509        VspBspTraversalData tData(mRoot,
510                                                          new PolygonContainer(polys),
[463]511                                                          0,
[482]512                                                          rays,
513                              ComputePvsSize(*rays),
[547]514                                                          prop,
[463]515                                                          geom);
[663]516
[1302]517        // evaluate the priority of this traversal data
[664]518        EvalPriority(tData);
[663]519
[578]520        // first node is kd node, i.e. an axis aligned box
[710]521        if (1)
[578]522        tData.mIsKdNode = true;
523        else
524                tData.mIsKdNode = false;
[562]525
[600]526        tQueue.push(tData);
[463]527
[609]528
[1563]529        mTotalCost = tData.mPvs * tData.mProbability / mBoundingBox.GetVolume();
[607]530        mTotalPvsSize = tData.mPvs;
531       
[1020]532        // first subdivison statistics
533        AddSubdivisionStats(1, 0, 0, mTotalCost, (float)mTotalPvsSize);
534   
[574]535        mBspStats.Start();
[667]536        cout << "Constructing vsp bsp tree ... \n";
[463]537
[1016]538        const long startTime = GetTime();       
539        // used for intermediate time measurements and progress
540        long interTime = GetTime();
541
[664]542        int nLeaves = 500;
543        int nViewCells = 500;
[587]544
[542]545        mOutOfMemory = false;
[612]546        mCreatedViewCells = 0;
547       
[600]548        while (!tQueue.empty())
[463]549        {
[600]550                tData = tQueue.top();
551            tQueue.pop();               
[463]552
[508]553                if (0 && !mOutOfMemory)
554                {
555                        float mem = GetMemUsage();
[478]556
[508]557                        if (mem > mMaxMemory)
558                        {
559                                mOutOfMemory = true;
560                                Debug << "memory limit reached: " << mem << endl;
561                        }
562                }
563
[587]564                // subdivide leaf node
[1016]565                const BspNode *r = Subdivide(tQueue, tData);
[463]566
567                if (r == mRoot)
[482]568                        Debug << "VSP BSP tree construction time spent at root: "
[542]569                                  << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
570
[654]571                if (mBspStats.Leaves() >= nLeaves)
[542]572                {
[612]573                        nLeaves += 500;
574
[574]575                        cout << "leaves=" << mBspStats.Leaves() << endl;
[542]576                        Debug << "needed "
[587]577                                  << TimeDiff(interTime, GetTime())*1e-3
[612]578                                  << " secs to create 500 view cells" << endl;
[542]579                        interTime = GetTime();
580                }
[612]581
[664]582                if (mCreatedViewCells >= nViewCells)
[612]583                {
584                        nViewCells += 500;
585
586                        cout << "generated " << mCreatedViewCells << " viewcells" << endl;
587                }
[463]588        }
589
[542]590        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
[1020]591        cout << "finished in " << TimeDiff(startTime, GetTime())*1e-3 << "secs" << endl;
[463]592
[574]593        mBspStats.Stop();
[463]594}
595
[508]596
[653]597
[654]598void VspBspTree::ConstructWithSplitQueue(const PolygonContainer &polys,
[653]599                                                                                          RayInfoContainer *rays)
600{
601        VspBspSplitQueue tQueue;
602
603        mRoot = new BspLeaf();
604
605        // constrruct root node geometry
606        BspNodeGeometry *geom = new BspNodeGeometry();
607        ConstructGeometry(mRoot, *geom);
608
609        const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();
610
611        VspBspTraversalData tData(mRoot,
612                                                          new PolygonContainer(polys),
613                                                          0,
614                                                          rays,
615                              ComputePvsSize(*rays),
616                                                          prop,
617                                                          geom);
618
[663]619
[1076]620        // first node is kd node, i.e. an axis aligned box
621        if (1)
622        tData.mIsKdNode = true;
623        else
624                tData.mIsKdNode = false;
625
[660]626        // compute first split candidate
[1233]627        VspBspSubdivisionCandidate splitCandidate;
[1145]628        splitCandidate.mParentData = tData;
[653]629
[1233]630        EvalSubdivisionCandidate(splitCandidate);
[1145]631
[653]632        tQueue.push(splitCandidate);
633
[1563]634        mTotalCost = tData.mPvs * tData.mProbability / mBoundingBox.GetVolume();
[653]635        mTotalPvsSize = tData.mPvs;
636       
[1020]637        // first subdivison statistics
638        AddSubdivisionStats(1, 0, 0, mTotalCost, (float)mTotalPvsSize);
639   
640    mBspStats.Start();
[667]641        cout << "Constructing vsp bsp tree ... \n";
[653]642
643        long startTime = GetTime();     
[666]644        int nLeaves = 500;
645        int nViewCells = 500;
[653]646
647        // used for intermediate time measurements and progress
648        long interTime = GetTime();     
649
650        mOutOfMemory = false;
651
652        mCreatedViewCells = 0;
653       
654        while (!tQueue.empty())
655        {
656                splitCandidate = tQueue.top();
657            tQueue.pop();               
658
[654]659                // cost ratio of cost decrease / totalCost
[1288]660                float costRatio = splitCandidate.mRenderCostDecr / mTotalCost;
[654]661
[655]662                //Debug << "cost ratio: " << costRatio << endl;
[654]663                if (costRatio < mTermMinGlobalCostRatio)
[1449]664                {
665                        ++ mBspStats.mGlobalCostMisses;
666                }
667
[653]668                if (0 && !mOutOfMemory)
669                {
670                        float mem = GetMemUsage();
671                        if (mem > mMaxMemory)
672                        {
673                                mOutOfMemory = true;
674                                Debug << "memory limit reached: " << mem << endl;
675                        }
676                }
677
678                // subdivide leaf node
679                BspNode *r = Subdivide(tQueue, splitCandidate);
680
681                if (r == mRoot)
[1449]682                {
[2069]683                        cout << "VSP BSP tree construction time spent at root: "
684                                 << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
[1449]685                }
[653]686
[654]687                if (mBspStats.Leaves() >= nLeaves)
[653]688                {
689                        nLeaves += 500;
690
691                        cout << "leaves=" << mBspStats.Leaves() << endl;
692                        Debug << "needed "
693                                  << TimeDiff(interTime, GetTime())*1e-3
694                                  << " secs to create 500 view cells" << endl;
695                        interTime = GetTime();
696                }
697
698                if (mCreatedViewCells == nViewCells)
699                {
700                        nViewCells += 500;
701                        cout << "generated " << mCreatedViewCells << " viewcells" << endl;
702                }
703        }
704
705        Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
706        cout << "finished\n";
707
708        mBspStats.Stop();
709}
710
711
[654]712bool VspBspTree::LocalTerminationCriteriaMet(const VspBspTraversalData &data) const
[463]713{
[482]714        return
[463]715                (((int)data.mRays->size() <= mTermMinRays) ||
[473]716                 (data.mPvs <= mTermMinPvs)   ||
[547]717                 (data.mProbability <= mTermMinProbability) ||
[535]718                 (data.GetAvgRayContribution() > mTermMaxRayContribution) ||
[463]719                 (data.mDepth >= mTermMaxDepth));
720}
721
[508]722
[1020]723void VspBspTree::AddSubdivisionStats(const int viewCells,
724                                                                         const float renderCostDecr,
725                                                                         const float splitCandidateCost,
726                                                                         const float totalRenderCost,
727                                                                         const float avgRenderCost)
728{
729        mSubdivisionStats
730                        << "#ViewCells\n" << viewCells << endl
731                        << "#RenderCostDecrease\n" << renderCostDecr << endl
[1233]732                        << "#SubdivisionCandidateCost\n" << splitCandidateCost << endl
[1020]733                        << "#TotalRenderCost\n" << totalRenderCost << endl
734                        << "#AvgRenderCost\n" << avgRenderCost << endl;
735}
736
737
[654]738bool VspBspTree::GlobalTerminationCriteriaMet(const VspBspTraversalData &data) const
739{
740        return
[1302]741                (0
742                || mOutOfMemory
[654]743                || (mBspStats.Leaves() >= mMaxViewCells)
[1449]744                || (mBspStats.mGlobalCostMisses >= mTermGlobalCostMissTolerance)
[654]745                 );
746}
747
748
[600]749BspNode *VspBspTree::Subdivide(VspBspTraversalQueue &tQueue,
[463]750                                                           VspBspTraversalData &tData)
751{
[473]752        BspNode *newNode = tData.mNode;
753
[654]754        if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
[473]755        {
756                PolygonContainer coincident;
[482]757
[473]758                VspBspTraversalData tFrontData;
759                VspBspTraversalData tBackData;
[612]760
[473]761                // create new interior node and two leaf nodes
762                // or return leaf as it is (if maxCostRatio missed)
[653]763                int splitAxis;
764                bool splitFurther = true;
765                int maxCostMisses = tData.mMaxCostMisses;
766               
767                Plane3 splitPlane;
[2017]768                BspLeaf *leaf = static_cast<BspLeaf *>(tData.mNode);
[726]769               
770                // choose next split plane
[653]771                if (!SelectPlane(splitPlane, leaf, tData, tFrontData, tBackData, splitAxis))
[473]772                {
[653]773                        ++ maxCostMisses;
774
775                        if (maxCostMisses > mTermMissTolerance)
776                        {
777                                // terminate branch because of max cost
778                                ++ mBspStats.maxCostNodes;
779                                splitFurther = false;
780                        }
781                }
782       
[726]783                // if this a valid split => subdivide this node further
[1193]784
785                if (splitFurther)
[653]786                {
787                        newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
788
789                        if (splitAxis < 3)
790                                ++ mBspStats.splits[splitAxis];
791                        else
792                                ++ mBspStats.polySplits;
793
[822]794                        // if it was a kd node (i.e., a box) and the split axis is axis aligned, it is still a kd node
[710]795                        tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && (splitAxis < 3));
[1076]796                       
[726]797                        tFrontData.mAxis = tBackData.mAxis = splitAxis;
[663]798
[653]799                        // how often was max cost ratio missed in this branch?
800                        tFrontData.mMaxCostMisses = maxCostMisses;
801                        tBackData.mMaxCostMisses = maxCostMisses;
802
[664]803                        EvalPriority(tFrontData);
804                        EvalPriority(tBackData);
[663]805
[726]806                        // evaluate subdivision stats
[600]807                        if (1)
[1145]808                                EvalSubdivisionStats(tData, tFrontData, tBackData);
809                       
[605]810
[473]811                        // push the children on the stack
[600]812                        tQueue.push(tFrontData);
813                        tQueue.push(tBackData);
[473]814
815                        // delete old leaf node
[482]816                        DEL_PTR(tData.mNode);
[473]817                }
818        }
[482]819
[478]820        //-- terminate traversal and create new view cell
[473]821        if (newNode->IsLeaf())
[463]822        {
[2017]823                BspLeaf *leaf = static_cast<BspLeaf *>(newNode);
[1020]824               
[547]825                BspViewCell *viewCell = new BspViewCell();
[463]826                leaf->SetViewCell(viewCell);
[487]827       
[2066]828                if (STORE_PVS)
829                {
830                        //////////
831                        //-- update pvs
[487]832
[2066]833                        int conSamp = 0;
834                        float sampCon = 0.0f;
[752]835
[2066]836                        AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
[487]837
[2066]838                        // update scalar pvs size lookup
839                        ObjectPvs &pvs = viewCell->GetPvs();
840                        mViewCellsManager->UpdateScalarPvsSize(viewCell, pvs.EvalPvsCost(), pvs.GetSize());
841               
842                        mBspStats.contributingSamples += conSamp;
843                        mBspStats.sampleContributions += (int)sampCon;
844                }
845
846                //////////
[487]847                //-- store additional info
[2066]848
[478]849                if (mStoreRays)
850                {
851                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
852                        for (it = tData.mRays->begin(); it != it_end; ++ it)
[639]853                        {
854                                (*it).mRay->Ref();                     
[478]855                                leaf->mVssRays.push_back((*it).mRay);
[639]856                        }
[478]857                }
[612]858
[564]859                // should I check here?
[1027]860                if (0 && !mViewCellsManager->CheckValidity(viewCell, 0,
861                        mViewCellsManager->GetMaxPvsSize()))
[547]862                {
863                        viewCell->SetValid(false);
864                        leaf->SetTreeValid(false);
865                        PropagateUpValidity(leaf);
[463]866
[574]867                        ++ mBspStats.invalidLeaves;
[547]868                }
869               
[1551]870                viewCell->mLeaves.push_back(leaf);
[547]871
872                if (mUseAreaForPvs)
873                        viewCell->SetArea(tData.mProbability);
874                else
875                        viewCell->SetVolume(tData.mProbability);
876
877                leaf->mProbability = tData.mProbability;
878
[1020]879                // finally evaluate stats until this leaf
[1027]880                if (0)
881                        EvaluateLeafStats(tData);               
[463]882        }
[482]883
[473]884        //-- cleanup
[478]885        tData.Clear();
[463]886
[472]887        return newNode;
[463]888}
889
[1027]890
891// subdivide node using a split plane queue
[653]892BspNode *VspBspTree::Subdivide(VspBspSplitQueue &tQueue,
[1233]893                                                           VspBspSubdivisionCandidate &splitCandidate)
[653]894{
895        VspBspTraversalData &tData = splitCandidate.mParentData;
896
897        BspNode *newNode = tData.mNode;
898
[654]899        if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
[666]900        {       
[653]901                PolygonContainer coincident;
902
903                VspBspTraversalData tFrontData;
904                VspBspTraversalData tBackData;
905
[1586]906                ////////////////////
[653]907                //-- continue subdivision
[726]908               
[653]909                // create new interior node and two leaf node
910                const Plane3 splitPlane = splitCandidate.mSplitPlane;
[654]911                               
[653]912                newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
[666]913       
[660]914                const int splitAxis = splitCandidate.mSplitAxis;
915                const int maxCostMisses = splitCandidate.mMaxCostMisses;
916
[654]917                if (splitAxis < 3)
918                        ++ mBspStats.splits[splitAxis];
919                else
920                        ++ mBspStats.polySplits;
[653]921
[710]922                tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && (splitAxis < 3));
[726]923                tFrontData.mAxis = tBackData.mAxis = splitAxis;
924
[654]925                // how often was max cost ratio missed in this branch?
926                tFrontData.mMaxCostMisses = maxCostMisses;
927                tBackData.mMaxCostMisses = maxCostMisses;
[664]928                       
[1020]929                // statistics
[653]930                if (1)
931                {
[2070]932                        const float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
933                        const float cBack = (float)tBackData.mPvs * tBackData.mProbability;
934                        const float cData = (float)tData.mPvs * tData.mProbability;
[675]935                       
[2070]936                        const float costDecr =
[1563]937                                (cFront + cBack - cData) / mBoundingBox.GetVolume();
[653]938
939                        mTotalCost += costDecr;
940                        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
941
[1020]942                        AddSubdivisionStats(mBspStats.Leaves(),
943                                                                -costDecr, 
944                                                                splitCandidate.GetPriority(),
945                                                                mTotalCost,
946                                                                (float)mTotalPvsSize / (float)mBspStats.Leaves());
[653]947                }
948
[2070]949                ////////////
[653]950                //-- push the new split candidates on the stack
[2070]951
[1233]952                VspBspSubdivisionCandidate frontCandidate;
[1145]953                frontCandidate.mParentData = tFrontData;
954
[1233]955                VspBspSubdivisionCandidate backCandidate;
[1145]956                backCandidate.mParentData = tBackData;
[653]957
[1233]958                EvalSubdivisionCandidate(frontCandidate);
959                EvalSubdivisionCandidate(backCandidate);
[654]960       
[2071]961                cout << "f cost: " << frontCandidate.mPriority << " " << frontCandidate.mRenderCostDecr << endl;
962                cout << "b cost: " << backCandidate.mPriority << " " << backCandidate.mRenderCostDecr << endl;
[653]963                tQueue.push(frontCandidate);
964                tQueue.push(backCandidate);
[666]965       
[653]966                // delete old leaf node
967                DEL_PTR(tData.mNode);
968        }
969
[654]970
[1551]971        //////////////////
[653]972        //-- terminate traversal and create new view cell
[1632]973
[653]974        if (newNode->IsLeaf())
975        {
[2017]976                BspLeaf *leaf = static_cast<BspLeaf *>(newNode);
[1020]977
[653]978                BspViewCell *viewCell = new BspViewCell();
[710]979        leaf->SetViewCell(viewCell);
[653]980               
[2066]981                if (STORE_PVS)
982                {
983                        /////////
984                        //-- update pvs
[653]985
[2066]986            int conSamp = 0;
987                        float sampCon = 0.0f;
[752]988
[2066]989                        AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
990               
991                        // update scalar pvs size value
992                        ObjectPvs &pvs = viewCell->GetPvs();
993                        mViewCellsManager->UpdateScalarPvsSize(viewCell,
994                                                                                                   pvs.EvalPvsCost(),
995                                                                                                   pvs.GetSize());
996                       
997                        mBspStats.contributingSamples += conSamp;
998                        mBspStats.sampleContributions += (int)sampCon;
999                }
[653]1000
[1551]1001                viewCell->mLeaves.push_back(leaf);
1002
1003                ///////////
[653]1004                //-- store additional info
[1632]1005
[653]1006                if (mStoreRays)
1007                {
1008                        RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
1009                        for (it = tData.mRays->begin(); it != it_end; ++ it)
1010                        {
1011                                (*it).mRay->Ref();                     
1012                                leaf->mVssRays.push_back((*it).mRay);
1013                        }
1014                }
[1551]1015       
[653]1016                if (mUseAreaForPvs)
1017                        viewCell->SetArea(tData.mProbability);
1018                else
1019                        viewCell->SetVolume(tData.mProbability);
1020
[675]1021        leaf->mProbability = tData.mProbability;
[653]1022
[2070]1023                // finally evaluate stats for this leaf
1024                if (1)
[1027]1025                        EvaluateLeafStats(tData);               
[653]1026        }
1027
1028        //-- cleanup
1029        tData.Clear();
1030
1031        return newNode;
1032}
1033
1034
[664]1035void VspBspTree::EvalPriority(VspBspTraversalData &tData) const
1036{
[735]1037    switch (mNodePriorityQueueType)
1038        {
1039        case BREATH_FIRST:
[734]1040                tData.mPriority = (float)-tData.mDepth;
[735]1041                break;
1042        case DEPTH_FIRST:
[734]1043                tData.mPriority = (float)tData.mDepth;
[735]1044                break;
1045        default:
[734]1046                tData.mPriority = tData.mPvs * tData.mProbability;
[735]1047                //Debug << "priority: " << tData.mPriority << endl;
1048                break;
1049        }
[664]1050}
1051
1052
[1233]1053void VspBspTree::EvalSubdivisionCandidate(VspBspSubdivisionCandidate &splitCandidate)
[1145]1054{
1055        VspBspTraversalData frontData;
1056        VspBspTraversalData backData;
1057       
[2017]1058        BspLeaf *leaf = static_cast<BspLeaf *>(splitCandidate.mParentData.mNode);
[1145]1059       
1060        // compute locally best split plane
[1288]1061    const bool costRatioViolated =
1062                SelectPlane(splitCandidate.mSplitPlane,
1063                                        leaf,
1064                                        splitCandidate.mParentData,
1065                                        frontData,
1066                                        backData,
1067                                        splitCandidate.mSplitAxis);
[652]1068
[1288]1069        // max cost threshold violated?
1070        splitCandidate.mMaxCostMisses = costRatioViolated ?
1071                splitCandidate.mParentData.mMaxCostMisses :
1072                splitCandidate.mParentData.mMaxCostMisses + 1;
1073
[1145]1074        float oldRenderCost;
1075
1076        // compute global decrease in render cost
1077        const float renderCostDecr = EvalRenderCostDecrease(splitCandidate.mSplitPlane,
1078                                                                                                                splitCandidate.mParentData,
1079                                                                                                                oldRenderCost);
1080
1081        splitCandidate.mRenderCostDecr = renderCostDecr;
1082
1083        // TODO: geometry could be reused
1084        delete frontData.mGeometry;
1085        delete backData.mGeometry;
1086
1087        // set priority for queue
1088#if 0
1089        const float priority = (float)-data.mDepth;
1090#else   
1091
1092        // take render cost of node into account
1093        // otherwise danger of being stuck in a local minimum!!
1094        const float factor = mRenderCostDecreaseWeight;
1095        const float priority = factor * renderCostDecr + (1.0f - factor) * oldRenderCost;
1096#endif
1097       
1098        splitCandidate.mPriority = priority;
1099}
1100
1101
1102void VspBspTree::EvalSubdivisionStats(const VspBspTraversalData &tData,
1103                                                                          const VspBspTraversalData &tFrontData,
1104                                                                          const VspBspTraversalData &tBackData)
1105{
1106        const float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
1107        const float cBack = (float)tBackData.mPvs * tBackData.mProbability;
1108        const float cData = (float)tData.mPvs * tData.mProbability;
1109       
1110        const float costDecr =
[1563]1111                (cFront + cBack - cData) / mBoundingBox.GetVolume();
[1145]1112
1113        mTotalCost += costDecr;
1114        mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
1115
1116        AddSubdivisionStats(mBspStats.Leaves(),
1117                                                -costDecr,
1118                                                0,
1119                                                mTotalCost,
1120                                                (float)mTotalPvsSize / (float)mBspStats.Leaves());
1121}
1122
1123
[653]1124BspInterior *VspBspTree::SubdivideNode(const Plane3 &splitPlane,
1125                                                                           VspBspTraversalData &tData,
1126                                                                           VspBspTraversalData &frontData,
1127                                                                           VspBspTraversalData &backData,
1128                                                                           PolygonContainer &coincident)
[463]1129{
[2017]1130        BspLeaf *leaf = static_cast<BspLeaf *>(tData.mNode);
[508]1131       
[473]1132        //-- the front and back traversal data is filled with the new values
1133        frontData.mDepth = tData.mDepth + 1;
[508]1134        frontData.mPolygons = new PolygonContainer();
[473]1135        frontData.mRays = new RayInfoContainer();
[508]1136       
[473]1137        backData.mDepth = tData.mDepth + 1;
[508]1138        backData.mPolygons = new PolygonContainer();
[473]1139        backData.mRays = new RayInfoContainer();
[508]1140       
[653]1141
1142        //-- subdivide rays
[652]1143        SplitRays(splitPlane,
[482]1144                          *tData.mRays,
1145                          *frontData.mRays,
[463]1146                          *backData.mRays);
[482]1147
[463]1148
[472]1149        // compute pvs
[463]1150        frontData.mPvs = ComputePvsSize(*frontData.mRays);
1151        backData.mPvs = ComputePvsSize(*backData.mRays);
1152
[508]1153        // split front and back node geometry and compute area
[547]1154       
1155        // if geometry was not already computed
[602]1156        if (!frontData.mGeometry && !backData.mGeometry)
[463]1157        {
[547]1158                frontData.mGeometry = new BspNodeGeometry();
1159                backData.mGeometry = new BspNodeGeometry();
[482]1160
[547]1161                tData.mGeometry->SplitGeometry(*frontData.mGeometry,
1162                                                                           *backData.mGeometry,
[652]1163                                                                           splitPlane,
[1563]1164                                                                           mBoundingBox,
[679]1165                                                                           //0.0f);
1166                                                                           mEpsilon);
[508]1167               
[547]1168                if (mUseAreaForPvs)
1169                {
1170                        frontData.mProbability = frontData.mGeometry->GetArea();
1171                        backData.mProbability = backData.mGeometry->GetArea();
[508]1172                }
[547]1173                else
1174                {
1175                        frontData.mProbability = frontData.mGeometry->GetVolume();
[654]1176                        backData.mProbability = tData.mProbability - frontData.mProbability;
[676]1177
[744]1178                        // should never come here: wrong volume !!!
[676]1179                        if (0)
1180                        {
[744]1181                                if (frontData.mProbability < -0.00001)
1182                                        Debug << "fatal error f: " << frontData.mProbability << endl;
1183                                if (backData.mProbability < -0.00001)
1184                                        Debug << "fatal error b: " << backData.mProbability << endl;
1185
1186                                // clamp because of precision issues
[676]1187                                if (frontData.mProbability < 0) frontData.mProbability = 0;
1188                                if (backData.mProbability < 0) backData.mProbability = 0;
1189                        }
[547]1190                }
[463]1191        }
[663]1192
[547]1193       
[652]1194    // subdivide polygons
1195        SplitPolygons(splitPlane,
1196                                  *tData.mPolygons,
1197                      *frontData.mPolygons,
1198                                  *backData.mPolygons,
1199                                  coincident);
[463]1200
[652]1201
1202
[653]1203        ///////////////////////////////////////
1204        // subdivide further
[652]1205
[711]1206        // store maximal and minimal depth
1207        if (tData.mDepth > mBspStats.maxDepth)
1208        {
1209                Debug << "max depth increases to " << tData.mDepth << " at " << mBspStats.Leaves() << " leaves" << endl;
1210                mBspStats.maxDepth = tData.mDepth;
1211        }
1212
[652]1213        mBspStats.nodes += 2;
1214
[711]1215   
[652]1216        BspInterior *interior = new BspInterior(splitPlane);
1217
[1715]1218#ifdef GTP_DEBUG
[652]1219        Debug << interior << endl;
1220#endif
1221
[711]1222
[463]1223        //-- create front and back leaf
1224
1225        BspInterior *parent = leaf->GetParent();
1226
1227        // replace a link from node's parent
[487]1228        if (parent)
[463]1229        {
1230                parent->ReplaceChildLink(leaf, interior);
1231                interior->SetParent(parent);
1232        }
1233        else // new root
1234        {
1235                mRoot = interior;
1236        }
1237
1238        // and setup child links
1239        interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
[482]1240
[463]1241        frontData.mNode = interior->GetFront();
1242        backData.mNode = interior->GetBack();
[473]1243
[650]1244        interior->mTimeStamp = mTimeStamp ++;
[652]1245       
[711]1246
[463]1247        //DEL_PTR(leaf);
1248        return interior;
1249}
1250
[508]1251
[463]1252void VspBspTree::AddToPvs(BspLeaf *leaf,
[482]1253                                                  const RayInfoContainer &rays,
[556]1254                                                  float &sampleContributions,
[463]1255                                                  int &contributingSamples)
1256{
[1002]1257        sampleContributions = 0;
1258        contributingSamples = 0;
[556]1259 
[1002]1260        RayInfoContainer::const_iterator it, it_end = rays.end();
[556]1261 
[1002]1262        ViewCellLeaf *vc = leaf->GetViewCell();
[556]1263 
[1002]1264        // add contributions from samples to the PVS
1265        for (it = rays.begin(); it != it_end; ++ it)
[463]1266        {
[1002]1267                float sc = 0.0f;
1268                VssRay *ray = (*it).mRay;
1269
1270                bool madeContrib = false;
1271                float contribution;
1272
1273                if (ray->mTerminationObject)
1274                {
[2071]1275                        //cout << "m";
[1002]1276                        if (vc->AddPvsSample(ray->mTerminationObject, ray->mPdf, contribution))
1277                                madeContrib = true;
1278                        sc += contribution;
1279                }
[556]1280         
[1006]1281                // only count termination objects?
[1020]1282                if (COUNT_ORIGIN_OBJECTS && ray->mOriginObject)
[1002]1283                {
1284                        if (vc->AddPvsSample(ray->mOriginObject, ray->mPdf, contribution))
1285                                madeContrib = true;
[1020]1286
[1002]1287                        sc += contribution;
1288                }
[556]1289               
[1047]1290                sampleContributions += sc;
1291               
1292                if (madeContrib)
1293                        ++ contributingSamples;
[463]1294        }
1295}
1296
[580]1297
[1233]1298void VspBspTree::SortSubdivisionCandidates(const RayInfoContainer &rays,
[710]1299                                                                         const int axis,
1300                                                                         float minBand,
1301                                                                         float maxBand)
[463]1302{
[1233]1303        mLocalSubdivisionCandidates->clear();
[463]1304
[2070]1305        const int requestedSize = 2 * (int)(rays.size());
1306
[480]1307        // creates a sorted split candidates array
[1233]1308        if (mLocalSubdivisionCandidates->capacity() > 500000 &&
1309                requestedSize < (int)(mLocalSubdivisionCandidates->capacity() / 10) )
[480]1310        {
[1233]1311        delete mLocalSubdivisionCandidates;
1312                mLocalSubdivisionCandidates = new vector<SortableEntry>;
[480]1313        }
[463]1314
[1233]1315        mLocalSubdivisionCandidates->reserve(requestedSize);
[480]1316
[710]1317        if (0)
[1137]1318        {       // float values => don't compare with exact values
[710]1319                minBand += Limits::Small;
1320                maxBand -= Limits::Small;
1321        }
1322
[480]1323        // insert all queries
[612]1324        for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
[473]1325        {
[612]1326                const bool positive = (*ri).mRay->HasPosDir(axis);
[1137]1327                float pos = (*ri).ExtrapOrigin(axis);
1328
[710]1329                // clamp to min / max band
1330                if (0) ClipValue(pos, minBand, maxBand);
1331               
[2070]1332                mLocalSubdivisionCandidates->
1333                        push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
[810]1334                                                                        pos, (*ri).mRay));
[480]1335
[710]1336                pos = (*ri).ExtrapTermination(axis);
[1137]1337
[710]1338                // clamp to min / max band
1339                if (0) ClipValue(pos, minBand, maxBand);
1340
[2070]1341                mLocalSubdivisionCandidates->
1342                        push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
[810]1343                                                                        pos, (*ri).mRay));
[473]1344        }
[480]1345
[2070]1346        sort(mLocalSubdivisionCandidates->begin(), mLocalSubdivisionCandidates->end());
[463]1347}
1348
[580]1349
[480]1350float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
1351                                                                                  const AxisAlignedBox3 &box,
1352                                                                                  const int pvsSize,
[710]1353                                                                                  const int axis,
[480]1354                                          float &position)
[463]1355{
[1149]1356        RayInfoContainer usedRays;
[1147]1357
[1727]1358        if (mMaxTests < (int)rays.size())
[1147]1359        {
[1149]1360                GetRayInfoSets(rays, mMaxTests, usedRays);
[1147]1361        }
1362        else
[1149]1363        {
[1147]1364                usedRays = rays;
[1149]1365        }
[1147]1366
[710]1367        const float minBox = box.Min(axis);
1368        const float maxBox = box.Max(axis);
[822]1369
[710]1370        const float sizeBox = maxBox - minBox;
[480]1371
[801]1372        const float minBand = minBox + mMinBand * sizeBox;
1373        const float maxBand = minBox + mMaxBand * sizeBox;
[710]1374
[1233]1375        SortSubdivisionCandidates(usedRays, axis, minBand, maxBand);
[710]1376
[1580]1377        //////////////////
[463]1378        // go through the lists, count the number of objects left and right
1379        // and evaluate the following cost funcion:
[480]1380        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
1381
[1789]1382        float pvsl = 0;
1383        float pvsr = (float)pvsSize;
[480]1384
[1789]1385        float pvsBack = pvsl;
1386        float pvsFront = pvsr;
[612]1387
1388        float sum = (float)pvsSize * sizeBox;
[463]1389        float minSum = 1e20f;
[822]1390       
[710]1391        // if no border can be found, take mid split
1392        position = minBox + 0.5f * sizeBox;
[822]1393       
1394        // the relative cost ratio
[1074]1395        float ratio = 99999999.0f;
[822]1396        bool splitPlaneFound = false;
[710]1397
[480]1398        Intersectable::NewMail();
[1147]1399        RayInfoContainer::const_iterator ri, ri_end = usedRays.end();
[612]1400
1401        // set all object as belonging to the front pvs
[1149]1402        for(ri = usedRays.begin(); ri != ri_end; ++ ri)
[463]1403        {
[612]1404                Intersectable *oObject = (*ri).mRay->mOriginObject;
1405                Intersectable *tObject = (*ri).mRay->mTerminationObject;
1406
[1020]1407                if (COUNT_ORIGIN_OBJECTS && oObject)
[463]1408                {
[612]1409                        if (!oObject->Mailed())
[482]1410                        {
[612]1411                                oObject->Mail();
1412                                oObject->mCounter = 1;
[482]1413                        }
[612]1414                        else
1415                        {
1416                                ++ oObject->mCounter;
1417                        }
[463]1418                }
[1020]1419
[612]1420                if (tObject)
1421                {
1422                        if (!tObject->Mailed())
1423                        {
1424                                tObject->Mail();
1425                                tObject->mCounter = 1;
1426                        }
1427                        else
1428                        {
1429                                ++ tObject->mCounter;
1430                        }
1431                }
[480]1432        }
1433
1434        Intersectable::NewMail();
[1233]1435        vector<SortableEntry>::const_iterator ci, ci_end = mLocalSubdivisionCandidates->end();
[612]1436
[1233]1437        for (ci = mLocalSubdivisionCandidates->begin(); ci != ci_end; ++ ci)
[480]1438        {
1439                VssRay *ray;
[612]1440                ray = (*ci).ray;
1441               
1442                Intersectable *oObject = ray->mOriginObject;
1443                Intersectable *tObject = ray->mTerminationObject;
1444               
[480]1445                switch ((*ci).type)
[463]1446                {
[480]1447                        case SortableEntry::ERayMin:
1448                                {
[1020]1449                                        if (COUNT_ORIGIN_OBJECTS && oObject && !oObject->Mailed())
[480]1450                                        {
[612]1451                                                oObject->Mail();
1452                                                ++ pvsl;
[480]1453                                        }
[1020]1454
[612]1455                                        if (tObject && !tObject->Mailed())
1456                                        {
1457                                                tObject->Mail();
1458                                                ++ pvsl;
1459                                        }
[1020]1460
[480]1461                                        break;
1462                                }
1463                        case SortableEntry::ERayMax:
1464                                {
[1020]1465                                        if (COUNT_ORIGIN_OBJECTS && oObject)
[612]1466                                        {
1467                                                if (-- oObject->mCounter == 0)
1468                                                        -- pvsr;
1469                                        }
[463]1470
[612]1471                                        if (tObject)
[480]1472                                        {
[612]1473                                                if (-- tObject->mCounter == 0)
1474                                                        -- pvsr;
[480]1475                                        }
1476
1477                                        break;
1478                                }
1479                }
[822]1480               
1481               
[1020]1482                // Note: we compare size of bounding boxes of front and back side because
1483                // of efficiency reasons (otherwise a new geometry would have to be computed
1484                // in each step and incremential evaluation would be difficult.
1485                // but then errors happen if the geometry is not an axis aligned box
1486                // (i.e., if a geometry aligned split was taken before)
1487                // question: is it sufficient to make this approximation?
[710]1488                if (((*ci).value >= minBand) && ((*ci).value <= maxBand))
[480]1489                {
[612]1490                        sum = pvsl * ((*ci).value - minBox) + pvsr * (maxBox - (*ci).value);
[480]1491
[1147]1492                        float currentPos;
1493                       
1494                        // HACK: current positition is BETWEEN visibility events
[1149]1495                        if (0 && ((ci + 1) != ci_end))
1496                        {
[1147]1497                                currentPos = ((*ci).value + (*(ci + 1)).value) * 0.5f;
[1149]1498                        }
[1147]1499                        else
[1149]1500                currentPos = (*ci).value;                       
[1147]1501
[822]1502                        //Debug  << "pos=" << (*ci).value << "\t pvs=(" <<  pvsl << "," << pvsr << ")" << endl;
1503                        //Debug << "cost= " << sum << endl;
[480]1504
1505                        if (sum < minSum)
[463]1506                        {
[822]1507                                splitPlaneFound = true;
1508
[463]1509                                minSum = sum;
[1147]1510                                position = currentPos;
[710]1511                               
[612]1512                                pvsBack = pvsl;
1513                                pvsFront = pvsr;
[463]1514                        }
1515                }
1516        }
[710]1517       
[1767]1518        ///////
1519        //-- compute cost
1520
[1772]1521        const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
1522        const float upperPvsLimit = (float)mViewCellsManager->GetMaxPvsSize();
[612]1523
1524        const float pOverall = sizeBox;
1525
1526        const float pBack = position - minBox;
1527        const float pFront = maxBox - position;
1528       
[1789]1529        const float penaltyOld = EvalPvsPenalty((float)pvsSize, lowerPvsLimit, upperPvsLimit);
[612]1530    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1531        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1532       
1533        const float oldRenderCost = penaltyOld * pOverall;
1534        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1535
[822]1536        if (splitPlaneFound)
1537        {
1538                ratio = mPvsFactor * newRenderCost / (oldRenderCost + Limits::Small);
1539        }
1540        //if (axis != 1)
1541        //Debug << "axis=" << axis << " costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1542         //    <<"\t pb=(" << pvsBack << ")\t pf=(" << pvsFront << ")" << endl;
[612]1543
[480]1544        return ratio;
[463]1545}
1546
[480]1547
[482]1548float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
[491]1549                                                                                 const VspBspTraversalData &tData,
[495]1550                                                                                 int &axis,
[508]1551                                                                                 BspNodeGeometry **frontGeom,
1552                                                                                 BspNodeGeometry **backGeom,
[547]1553                                                                                 float &pFront,
1554                                                                                 float &pBack,
[710]1555                                                                                 const bool isKdNode)
[463]1556{
[508]1557        float nPosition[3];
1558        float nCostRatio[3];
[547]1559        float nProbFront[3];
1560        float nProbBack[3];
[508]1561
1562        BspNodeGeometry *nFrontGeom[3];
1563        BspNodeGeometry *nBackGeom[3];
1564
[822]1565        // set to NULL, so I can find out which gemetry was stored
[612]1566        for (int i = 0; i < 3; ++ i)
1567        {
1568                nFrontGeom[i] = NULL;
1569                nBackGeom[i] = NULL;
1570        }
1571
[545]1572        // create bounding box of node geometry
[463]1573        AxisAlignedBox3 box;
[710]1574               
[562]1575        //TODO: for kd split geometry already is box => only take minmax vertices
[551]1576        if (1)
[1147]1577        {       // get bounding box from geometry
[710]1578                tData.mGeometry->GetBoundingBox(box);
[509]1579        }
1580        else
1581        {
[710]1582                box.Initialize();
[509]1583                RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();
[480]1584
[509]1585                for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
1586                        box.Include((*ri).ExtrapTermination());
1587        }
[663]1588
[1147]1589
[663]1590        int sAxis = 0;
[822]1591        int bestAxis;
1592
[837]1593        // if max cost ratio is exceeded, take split along longest axis instead
1594        const float maxCostRatioForArbitraryAxis = 0.9f;
1595
[1006]1596        if (mUseDrivingAxisIfMaxCostViolated)
[822]1597                bestAxis = box.Size().DrivingAxis();
1598        else
1599                bestAxis = -1;
1600
[978]1601#if 0
1602        // maximum cost ratio for axis to be valid:
[822]1603        // if exceeded, spatial mid split is used instead
[978]1604        const maxCostRatioForHeur = 0.99f;
1605#endif
[822]1606
[978]1607        // if we use some kind of specialised fixed axis
[837]1608    const bool useSpecialAxis =
1609                mOnlyDrivingAxis || mUseRandomAxis || mCirculatingAxis;
[663]1610
[801]1611        if (mUseRandomAxis)
[663]1612                sAxis = Random(3);
[801]1613        else if (mCirculatingAxis)
[726]1614                sAxis = (tData.mAxis + 1) % 3;
[822]1615        else if (mOnlyDrivingAxis)
[801]1616                sAxis = box.Size().DrivingAxis();
1617
[664]1618               
[670]1619        //Debug << "use special axis: " << useSpecialAxis << endl;
1620        //Debug << "axis: " << sAxis << " drivingaxis: " << box.Size().DrivingAxis();
[1076]1621       
[822]1622        for (axis = 0; axis < 3 ; ++ axis)
[480]1623        {
[663]1624                if (!useSpecialAxis || (axis == sAxis))
[463]1625                {
[822]1626                        if (mUseCostHeuristics)
[480]1627                        {
[1106]1628                                //-- place split plane using heuristics
[822]1629                                nCostRatio[axis] =
1630                                        BestCostRatioHeuristics(*tData.mRays,
1631                                                                                    box,
1632                                                                                        tData.mPvs,
1633                                                                                        axis,
1634                                                                                        nPosition[axis]);                       
1635                        }
[1106]1636                        else
[1767]1637                        {       
1638                                //-- split plane position is spatial median
[480]1639                                nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
[543]1640                                Vector3 normal(0,0,0); normal[axis] = 1.0f;
[1076]1641                               
[547]1642                                // allows faster split because we have axis aligned kd tree boxes
[822]1643                                if (isKdNode)
[542]1644                                {
1645                                        nCostRatio[axis] = EvalAxisAlignedSplitCost(tData,
1646                                                                                                                                box,
1647                                                                                                                                axis,
[547]1648                                                                                                                                nPosition[axis],
1649                                                                                                                                nProbFront[axis],
1650                                                                                                                                nProbBack[axis]);
[542]1651                                       
[710]1652                                        // create back geometry from box
[1767]1653
[1047]1654                                        // NOTE: the geometry is returned from the function so we
1655                                        // don't have to recompute it when possible
[1767]1656                                        Vector3 pos;
1657                                       
[543]1658                                        pos = box.Max(); pos[axis] = nPosition[axis];
[547]1659                                        AxisAlignedBox3 bBox(box.Min(), pos);
[1767]1660                                       
[710]1661                                        PolygonContainer fPolys;
1662                                        bBox.ExtractPolys(fPolys);
1663
1664                                        nBackGeom[axis] = new BspNodeGeometry(fPolys);
1665       
[1767]1666                                        ////////////
[822]1667                                        //-- create front geometry from box
[1767]1668
[542]1669                                        pos = box.Min(); pos[axis] = nPosition[axis];
[547]1670                                        AxisAlignedBox3 fBox(pos, box.Max());
[710]1671
1672                                        PolygonContainer bPolys;
1673                                        fBox.ExtractPolys(bPolys);
1674                                        nFrontGeom[axis] = new BspNodeGeometry(bPolys);
[542]1675                                }
1676                                else
1677                                {
[710]1678                                        nFrontGeom[axis] = new BspNodeGeometry();
1679                                        nBackGeom[axis] = new BspNodeGeometry();
1680
[573]1681                                        nCostRatio[axis] =
1682                                                EvalSplitPlaneCost(Plane3(normal, nPosition[axis]),
1683                                                                                   tData, *nFrontGeom[axis], *nBackGeom[axis],
1684                                                                                   nProbFront[axis], nProbBack[axis]);
[542]1685                                }
[480]1686                        }
[822]1687                                               
[978]1688                       
[1006]1689                        if (mUseDrivingAxisIfMaxCostViolated)
[480]1690                        {
[978]1691                                // we take longest axis split if cost ratio exceeds threshold
1692                                if (nCostRatio[axis] < min(maxCostRatioForArbitraryAxis, nCostRatio[bestAxis]))
[822]1693                                {
1694                                        bestAxis = axis;
1695                                }
[1072]1696                                /*else if (nCostRatio[axis] < nCostRatio[bestAxis])
1697                                {
[978]1698                                        Debug << "taking split along longest axis (" << bestAxis << ") instead of  (" << axis << ")" << endl;
[1072]1699                                }*/
1700
[978]1701                        }
1702                        else
1703                        {
1704                                if (bestAxis == -1)
[822]1705                                {
1706                                        bestAxis = axis;
1707                                }
[978]1708                                else if (nCostRatio[axis] < nCostRatio[bestAxis])
[822]1709                                {
1710                                        bestAxis = axis;
1711                                }
[978]1712                        }
[463]1713                }
1714        }
1715
[1767]1716        //////////
[495]1717        //-- assign values
[1072]1718
[495]1719        axis = bestAxis;
[547]1720        pFront = nProbFront[bestAxis];
1721        pBack = nProbBack[bestAxis];
[495]1722
[508]1723        // assign best split nodes geometry
1724        *frontGeom = nFrontGeom[bestAxis];
1725        *backGeom = nBackGeom[bestAxis];
[542]1726
[508]1727        // and delete other geometry
[612]1728        DEL_PTR(nFrontGeom[(bestAxis + 1) % 3]);
1729        DEL_PTR(nBackGeom[(bestAxis + 2) % 3]);
[495]1730
[508]1731        //-- split plane
1732    Vector3 normal(0,0,0); normal[bestAxis] = 1;
[480]1733        plane = Plane3(normal, nPosition[bestAxis]);
[508]1734
[822]1735        //Debug << "best axis: " << bestAxis << " pos " << nPosition[bestAxis] << endl;
[1072]1736
[480]1737        return nCostRatio[bestAxis];
[463]1738}
1739
[480]1740
[508]1741bool VspBspTree::SelectPlane(Plane3 &bestPlane,
1742                                                         BspLeaf *leaf,
[652]1743                                                         VspBspTraversalData &data,                                                     
[508]1744                                                         VspBspTraversalData &frontData,
[612]1745                                                         VspBspTraversalData &backData,
1746                                                         int &splitAxis)
[491]1747{
[810]1748        // HACK matt: subdivide regularily to certain depth
[1072]1749        if (data.mDepth < 0)    // question matt: why depth < 0 ?
[801]1750        {
[1072]1751                cout << "depth: " << data.mDepth << endl;
1752
[801]1753                // return axis aligned split
1754                AxisAlignedBox3 box;
1755                box.Initialize();
1756       
1757                // create bounding box of region
1758                data.mGeometry->GetBoundingBox(box);
1759       
1760                const int axis = box.Size().DrivingAxis();
1761                const Vector3 position = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1762
1763                Vector3 norm(0,0,0); norm[axis] = 1.0f;
1764                bestPlane = Plane3(norm, position);
1765                splitAxis = axis;
[1072]1766
[801]1767                return true;
1768        }
1769
[508]1770        // simplest strategy: just take next polygon
1771        if (mSplitPlaneStrategy & RANDOM_POLYGON)
[491]1772        {
[508]1773        if (!data.mPolygons->empty())
1774                {
1775                        const int randIdx =
1776                                (int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
1777                        Polygon3 *nextPoly = (*data.mPolygons)[randIdx];
[491]1778
[508]1779                        bestPlane = nextPoly->GetSupportingPlane();
1780                        return true;
1781                }
[491]1782        }
1783
[508]1784        //-- use heuristics to find appropriate plane
[491]1785
[508]1786        // intermediate plane
1787        Plane3 plane;
1788        float lowestCost = MAX_FLOAT;
[517]1789       
1790        // decides if the first few splits should be only axisAligned
[508]1791        const bool onlyAxisAligned  =
1792                (mSplitPlaneStrategy & AXIS_ALIGNED) &&
1793                ((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
1794                ((int)data.GetAvgRayContribution() < mTermMaxRayContriForAxisAligned);
1795       
1796        const int limit = onlyAxisAligned ? 0 :
1797                Min((int)data.mPolygons->size(), mMaxPolyCandidates);
[491]1798
[508]1799        float candidateCost;
[491]1800
[508]1801        int maxIdx = (int)data.mPolygons->size();
[491]1802
[508]1803        for (int i = 0; i < limit; ++ i)
[491]1804        {
[562]1805                // the already taken candidates are stored behind maxIdx
1806                // => assure that no index is taken twice
[508]1807                const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
1808                Polygon3 *poly = (*data.mPolygons)[candidateIdx];
[491]1809
[508]1810                // swap candidate to the end to avoid testing same plane
1811                std::swap((*data.mPolygons)[maxIdx], (*data.mPolygons)[candidateIdx]);
1812                //Polygon3 *poly = (*data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
[491]1813
[508]1814                // evaluate current candidate
1815                BspNodeGeometry fGeom, bGeom;
1816                float fArea, bArea;
1817                plane = poly->GetSupportingPlane();
[573]1818                candidateCost = EvalSplitPlaneCost(plane, data, fGeom, bGeom, fArea, bArea);
[491]1819               
[508]1820                if (candidateCost < lowestCost)
[491]1821                {
[508]1822                        bestPlane = plane;
1823                        lowestCost = candidateCost;
[491]1824                }
1825        }
1826
[1006]1827
[508]1828        //-- evaluate axis aligned splits
[1006]1829       
[508]1830        int axis;
1831        BspNodeGeometry *fGeom, *bGeom;
[547]1832        float pFront, pBack;
[491]1833
[653]1834        candidateCost = 99999999.0f;
[491]1835
[1006]1836        // as a variant, we take axis aligned split only if there is
1837        // more polygon available to guide the split
[653]1838        if (!mUsePolygonSplitIfAvailable || data.mPolygons->empty())
1839        {
1840                candidateCost = SelectAxisAlignedPlane(plane,
1841                                                                                           data,
1842                                                                                           axis,
1843                                                                                           &fGeom,
1844                                                                                           &bGeom,
1845                                                                                           pFront,
1846                                                                                           pBack,
1847                                                                                           data.mIsKdNode);     
1848        }
1849
[612]1850        splitAxis = 3;
[562]1851
[508]1852        if (candidateCost < lowestCost)
1853        {       
1854                bestPlane = plane;
1855                lowestCost = candidateCost;
[612]1856                splitAxis = axis;
[653]1857       
[542]1858                // assign already computed values
1859                // we can do this because we always save the
[562]1860                // computed values from the axis aligned splits         
[653]1861
[612]1862                if (fGeom && bGeom)
1863                {
1864                        frontData.mGeometry = fGeom;
1865                        backData.mGeometry = bGeom;
[547]1866       
[612]1867                        frontData.mProbability = pFront;
1868                        backData.mProbability = pBack;
1869                }
[508]1870        }
1871        else
[463]1872        {
[508]1873                DEL_PTR(fGeom);
1874                DEL_PTR(bGeom);
[463]1875        }
[678]1876   
[1715]1877#ifdef GTP_DEBUG
[508]1878        Debug << "plane lowest cost: " << lowestCost << endl;
[679]1879#endif
[508]1880
[801]1881        // exeeded relative max cost ratio
[508]1882        if (lowestCost > mTermMaxCostRatio)
[611]1883        {
[508]1884                return false;
[611]1885        }
[508]1886
1887        return true;
[463]1888}
1889
[480]1890
[473]1891Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
[482]1892{
[473]1893        const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
[482]1894
[473]1895        const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
1896        const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();
1897
1898        const Vector3 pt = (maxPt + minPt) * 0.5;
1899        const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());
1900
1901        return Plane3(normal, pt);
1902}
1903
[480]1904
[473]1905Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
[482]1906{
[473]1907        Vector3 pt[3];
[482]1908
[473]1909        int idx[3];
1910        int cmaxT = 0;
1911        int cminT = 0;
1912        bool chooseMin = false;
1913
1914        for (int j = 0; j < 3; ++ j)
1915        {
1916                idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));
[482]1917
[473]1918                if (idx[j] >= (int)rays.size())
1919                {
1920                        idx[j] -= (int)rays.size();
[482]1921
[473]1922                        chooseMin = (cminT < 2);
1923                }
1924                else
1925                        chooseMin = (cmaxT < 2);
1926
1927                RayInfo rayInf = rays[idx[j]];
1928                pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
[482]1929        }
[473]1930
1931        return Plane3(pt[0], pt[1], pt[2]);
1932}
1933
[580]1934
[473]1935Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
[482]1936{
[473]1937        Vector3 pt[3];
[482]1938
[473]1939        int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1940        int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1941
1942        // check if rays different
1943        if (idx1 == idx2)
1944                idx2 = (idx2 + 1) % (int)rays.size();
1945
1946        const RayInfo ray1 = rays[idx1];
1947        const RayInfo ray2 = rays[idx2];
1948
1949        // normal vector of the plane parallel to both lines
1950        const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));
1951
1952        // vector from line 1 to line 2
[479]1953        const Vector3 vd = ray2.ExtrapOrigin() - ray1.ExtrapOrigin();
[482]1954
[473]1955        // project vector on normal to get distance
1956        const float dist = DotProd(vd, norm);
1957
1958        // point on plane lies halfway between the two planes
1959        const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;
1960
1961        return Plane3(norm, planePt);
1962}
1963
[495]1964
[463]1965inline void VspBspTree::GenerateUniqueIdsForPvs()
1966{
[580]1967        Intersectable::NewMail(); sBackId = Intersectable::sMailId;
1968        Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
1969        Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
[463]1970}
1971
[495]1972
[652]1973float VspBspTree::EvalRenderCostDecrease(const Plane3 &candidatePlane,
[1145]1974                                                                                 const VspBspTraversalData &data,
1975                                                                                 float &normalizedOldRenderCost) const
[652]1976{
[729]1977        float pvsFront = 0;
1978        float pvsBack = 0;
1979        float totalPvs = 0;
[652]1980
1981        // probability that view point lies in back / front node
1982        float pOverall = data.mProbability;
1983        float pFront = 0;
1984        float pBack = 0;
1985
1986
1987        // create unique ids for pvs heuristics
1988        GenerateUniqueIdsForPvs();
1989       
1990        for (int i = 0; i < data.mRays->size(); ++ i)
1991        {
1992                RayInfo rayInf = (*data.mRays)[i];
1993
1994                float t;
1995                VssRay *ray = rayInf.mRay;
1996                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
1997
1998                // find front and back pvs for origing and termination object
1999                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
[1020]2000
[1577]2001                if (COUNT_ORIGIN_OBJECTS)
2002                {
[1020]2003                        AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
[1577]2004                }
[652]2005        }
2006
2007
2008        BspNodeGeometry geomFront;
2009        BspNodeGeometry geomBack;
2010
2011        // construct child geometry with regard to the candidate split plane
2012        data.mGeometry->SplitGeometry(geomFront,
2013                                                                  geomBack,
2014                                                                  candidatePlane,
[1563]2015                                                                  mBoundingBox,
[679]2016                                                                  //0.0f);
2017                                                                  mEpsilon);
[652]2018
2019        if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
2020        {
2021                pFront = geomFront.GetVolume();
2022                pBack = pOverall - pFront;
[675]2023
[729]2024                // something is wrong with the volume
[1027]2025                if (0 && ((pFront < 0.0) || (pBack < 0.0)))
[676]2026                {
[752]2027                        Debug << "ERROR in volume:\n"
2028                                  << "volume f :" << pFront << " b: " << pBack << " p: " << pOverall
2029                                  << ", real volume f: " << pFront << " b: " << geomBack.GetVolume()
2030                                  << ", #polygons f: " << geomFront.Size() << " b: " << geomBack.Size() << " p: " << data.mGeometry->Size() << endl;
[676]2031                }
[652]2032        }
2033        else
2034        {
2035                pFront = geomFront.GetArea();
2036                pBack = geomBack.GetArea();
2037        }
2038       
2039
2040        // -- pvs rendering heuristics
[1020]2041
2042        // upper and lower bounds
[1772]2043        const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
2044        const float upperPvsLimit = (float)mViewCellsManager->GetMaxPvsSize();
[652]2045
[1772]2046        const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
2047    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
2048        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
[652]2049                       
2050        const float oldRenderCost = pOverall * penaltyOld;
2051        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
2052
[1563]2053        const float renderCostDecrease = (oldRenderCost - newRenderCost) / mBoundingBox.GetVolume();
[752]2054       
[1006]2055        // take render cost of node into account to avoid being stuck in a local minimum
[1563]2056        normalizedOldRenderCost = oldRenderCost / mBoundingBox.GetVolume();
[1006]2057       
[1145]2058        return renderCostDecrease;
[652]2059}
2060
2061
[573]2062float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
2063                                                                         const VspBspTraversalData &data,
2064                                                                         BspNodeGeometry &geomFront,
2065                                                                         BspNodeGeometry &geomBack,
2066                                                                         float &pFront,
2067                                                                         float &pBack) const
[463]2068{
[729]2069        float totalPvs = 0;
2070        float pvsFront = 0;
2071        float pvsBack = 0;
[652]2072       
[1006]2073        // overall probability is used as normalizer
[463]2074        float pOverall = 0;
[1006]2075
2076        // probability that view point lies in back / front node
[547]2077        pFront = 0;
2078        pBack = 0;
[463]2079
[1020]2080        int numTests; // the number of tests
[508]2081
[1020]2082        // if random samples shold be taken instead of testing all the rays
2083        bool useRand;
[463]2084
2085        if ((int)data.mRays->size() > mMaxTests)
2086        {
2087                useRand = true;
[1020]2088                numTests = mMaxTests;
[463]2089        }
2090        else
2091        {
2092                useRand = false;
[1020]2093                numTests = (int)data.mRays->size();
[463]2094        }
[508]2095       
[1020]2096        // create unique ids for pvs heuristics
2097        GenerateUniqueIdsForPvs();
2098
2099        for (int i = 0; i < numTests; ++ i)
[463]2100        {
[508]2101                const int testIdx = useRand ?
2102                        (int)RandomValue(0, (Real)((int)data.mRays->size() - 1)) : i;
[463]2103                RayInfo rayInf = (*data.mRays)[testIdx];
2104
2105                float t;
[508]2106                VssRay *ray = rayInf.mRay;
[463]2107                const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
2108
[652]2109                // find front and back pvs for origing and termination object
2110                AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
[1020]2111
[1580]2112                if (COUNT_ORIGIN_OBJECTS)
2113                {
[1020]2114                        AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
[1580]2115                }
[652]2116        }
[508]2117
[652]2118        // construct child geometry with regard to the candidate split plane
[679]2119        bool splitSuccessFull = data.mGeometry->SplitGeometry(geomFront,
2120                                                                                                                  geomBack,
2121                                                                                                                  candidatePlane,
[1563]2122                                                                                                                  mBoundingBox,
[679]2123                                                                                                                  //0.0f);
2124                                                                                                                  mEpsilon);
[675]2125
[652]2126        pOverall = data.mProbability;
[463]2127
[652]2128        if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
2129        {
2130                pFront = geomFront.GetVolume();
2131                pBack = pOverall - pFront;
[675]2132               
[729]2133                // HACK: precision issues possible for unbalanced split => don't take this split!
[685]2134                if (1 &&
[682]2135                        (!splitSuccessFull || (pFront <= 0) || (pBack <= 0) ||
2136                        !geomFront.Valid() || !geomBack.Valid()))
[676]2137                {
[752]2138                        //Debug << "error f: " << pFront << " b: " << pBack << endl;
[1027]2139
2140                        // high penalty for degenerated / wrong split
[711]2141                        return 99999.9f;
[676]2142                }
[463]2143        }
[652]2144        else
[542]2145        {
[652]2146                pFront = geomFront.GetArea();
2147                pBack = geomBack.GetArea();
[542]2148        }
[652]2149       
[1772]2150        ////////
2151        //-- pvs rendering heuristics
[542]2152
[1773]2153        const float lowerPvsLimit = (float)mViewCellsManager->GetMinPvsSize();
2154        const float upperPvsLimit = (float)mViewCellsManager->GetMaxPvsSize();
[580]2155
[652]2156        // only render cost heuristics or combined with standard deviation
[1772]2157        const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
2158    const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
2159        const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
[579]2160                       
[652]2161        const float oldRenderCost = pOverall * penaltyOld;
2162        const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
[579]2163
[652]2164        float oldCost, newCost;
[580]2165
[652]2166        // only render cost
2167        if (1)
2168        {
2169                oldCost = oldRenderCost;
2170                newCost = newRenderCost;
2171        }
2172        else // also considering standard deviation
2173        {
2174                // standard deviation is difference of back and front pvs
2175                const float expectedCost = 0.5f * (penaltyFront + penaltyBack);
[580]2176
[652]2177                const float newDeviation = 0.5f *
2178                        fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);
[580]2179
[652]2180                const float oldDeviation = penaltyOld;
[580]2181
[652]2182                newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
2183                oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
2184        }
[580]2185
[1020]2186        const float cost = mPvsFactor * newCost / (oldCost + Limits::Small);
[607]2187               
[463]2188
[1715]2189#ifdef GTP_DEBUG
[474]2190        Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
[508]2191                  << " frontpvs: " << pvsFront << " pFront: " << pFront
2192                  << " backpvs: " << pvsBack << " pBack: " << pBack << endl << endl;
[600]2193        Debug << "cost: " << cost << endl;
[463]2194#endif
[482]2195
[652]2196        return cost;
[463]2197}
2198
[508]2199
[697]2200int VspBspTree::ComputeBoxIntersections(const AxisAlignedBox3 &box,
2201                                                                                ViewCellContainer &viewCells) const
2202{
2203        stack<bspNodePair> nodeStack;
2204        BspNodeGeometry *rgeom = new BspNodeGeometry();
2205
2206        ConstructGeometry(mRoot, *rgeom);
2207
2208        nodeStack.push(bspNodePair(mRoot, rgeom));
2209 
2210        ViewCell::NewMail();
2211
2212        while (!nodeStack.empty())
2213        {
2214                BspNode *node = nodeStack.top().first;
2215                BspNodeGeometry *geom = nodeStack.top().second;
2216                nodeStack.pop();
2217
2218                const int side = geom->ComputeIntersection(box);
2219               
2220                switch (side)
2221                {
2222                case -1:
2223                        // node geometry is contained in box
2224                        CollectViewCells(node, true, viewCells, true);
2225                        break;
2226
2227                case 0:
2228                        if (node->IsLeaf())
2229                        {
[2017]2230                                BspLeaf *leaf = static_cast<BspLeaf *>(node);
[697]2231                       
2232                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
2233                                {
2234                                        leaf->GetViewCell()->Mail();
2235                                        viewCells.push_back(leaf->GetViewCell());
2236                                }
2237                        }
2238                        else
2239                        {
[2017]2240                                BspInterior *interior = static_cast<BspInterior *>(node);
[697]2241                       
2242                                BspNode *first = interior->GetFront();
2243                                BspNode *second = interior->GetBack();
2244           
2245                                BspNodeGeometry *firstGeom = new BspNodeGeometry();
2246                                BspNodeGeometry *secondGeom = new BspNodeGeometry();
2247
2248                                geom->SplitGeometry(*firstGeom,
2249                                                                        *secondGeom,
2250                                                                        interior->GetPlane(),
[1563]2251                                                                        mBoundingBox,
[697]2252                                                                        //0.0000001f);
2253                                                                        mEpsilon);
2254
2255                                nodeStack.push(bspNodePair(first, firstGeom));
2256                                nodeStack.push(bspNodePair(second, secondGeom));
2257                        }
2258                       
2259                        break;
2260                default:
2261                        // default: cull
2262                        break;
2263                }
2264               
2265                DEL_PTR(geom);
2266               
2267        }
2268
2269        return (int)viewCells.size();
2270}
2271
2272
[580]2273float VspBspTree::EvalAxisAlignedSplitCost(const VspBspTraversalData &data,
2274                                                                                   const AxisAlignedBox3 &box,
2275                                                                                   const int axis,
2276                                                                                   const float &position,                                                                                 
2277                                                                                   float &pFront,
2278                                                                                   float &pBack) const
2279{
[729]2280        float pvsTotal = 0;
2281        float pvsFront = 0;
2282        float pvsBack = 0;
[580]2283       
2284        // create unique ids for pvs heuristics
2285        GenerateUniqueIdsForPvs();
2286
2287        const int pvsSize = data.mPvs;
[726]2288
[580]2289        RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
2290
2291        // this is the main ray classification loop!
2292        for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
2293        {
2294                // determine the side of this ray with respect to the plane
2295                float t;
2296                const int side = (*rit).ComputeRayIntersection(axis, position, t);
2297       
2298                AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
[1020]2299
[1765]2300                if (COUNT_ORIGIN_OBJECTS)
2301                {
[1020]2302                        AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
[1765]2303                }
[580]2304        }
2305
[1020]2306
[580]2307        //-- pvs heuristics
2308
[1020]2309        float pOverall = data.mProbability;
[580]2310
[1020]2311        // note: we use a simplified computation assuming that we always do a
2312        // spatial mid split   
2313       
[580]2314        if (!mUseAreaForPvs)
[1020]2315        {   
2316                // volume
[580]2317                pBack = pFront = pOverall * 0.5f;
2318#if 0
2319                // box length substitute for probability
2320                const float minBox = box.Min(axis);
2321                const float maxBox = box.Max(axis);
2322
2323                pBack = position - minBox;
2324                pFront = maxBox - position;
2325                pOverall = maxBox - minBox;
2326#endif
2327        }
2328        else //-- area substitute for probability
2329        {
2330                const int axis2 = (axis + 1) % 3;
2331                const int axis3 = (axis + 2) % 3;
2332
2333                const float faceArea =
2334                        (box.Max(axis2) - box.Min(axis2)) *
2335                        (box.Max(axis3) - box.Min(axis3));
2336
2337                pBack = pFront = pOverall * 0.5f + faceArea;
2338        }
2339
[1715]2340#ifdef GTP_DEBUG
[1302]2341        Debug << "axis: " << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
2342        Debug << "p: " << pFront << " " << pBack << " " << pOverall << endl;
[580]2343#endif
2344
2345       
2346        const float newCost = pvsBack * pBack + pvsFront * pFront;
2347        const float oldCost = (float)pvsSize * pOverall + Limits::Small;
2348
2349        return  (mCtDivCi + newCost) / oldCost;
2350}
2351
2352
[1076]2353inline void VspBspTree::AddObjToPvs(Intersectable *obj,
2354                                                                                 const int cf,
2355                                                                                 float &frontPvs,
2356                                                                                 float &backPvs,
2357                                                                                 float &totalPvs) const
[463]2358{
2359        if (!obj)
2360                return;
[1076]2361#if 0
[744]2362        const float renderCost = mViewCellsManager->EvalRenderCost(obj);
[1076]2363#else
2364        const int renderCost = 1;
2365#endif
[654]2366        // new object
[508]2367        if ((obj->mMailbox != sFrontId) &&
2368                (obj->mMailbox != sBackId) &&
2369                (obj->mMailbox != sFrontAndBackId))
2370        {
[744]2371                totalPvs += renderCost;
[508]2372        }
2373
[463]2374        // TODO: does this really belong to no pvs?
2375        //if (cf == Ray::COINCIDENT) return;
2376
2377        // object belongs to both PVS
2378        if (cf >= 0)
2379        {
[482]2380                if ((obj->mMailbox != sFrontId) &&
[463]2381                        (obj->mMailbox != sFrontAndBackId))
2382                {
[744]2383                        frontPvs += renderCost;
[508]2384               
[463]2385                        if (obj->mMailbox == sBackId)
[482]2386                                obj->mMailbox = sFrontAndBackId;
[463]2387                        else
[482]2388                                obj->mMailbox = sFrontId;
[463]2389                }
2390        }
[482]2391
[463]2392        if (cf <= 0)
2393        {
2394                if ((obj->mMailbox != sBackId) &&
2395                        (obj->mMailbox != sFrontAndBackId))
2396                {
[744]2397                        backPvs += renderCost;
[508]2398               
[463]2399                        if (obj->mMailbox == sFrontId)
[482]2400                                obj->mMailbox = sFrontAndBackId;
[463]2401                        else
[482]2402                                obj->mMailbox = sBackId;
[463]2403                }
2404        }
2405}
2406
[491]2407
[503]2408void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves,
2409                                                           const bool onlyUnmailed,
2410                                                           const int maxPvsSize) const
[463]2411{
2412        stack<BspNode *> nodeStack;
2413        nodeStack.push(mRoot);
[482]2414
2415        while (!nodeStack.empty())
[463]2416        {
2417                BspNode *node = nodeStack.top();
2418                nodeStack.pop();
[489]2419               
[482]2420                if (node->IsLeaf())
[463]2421                {
[490]2422                        // test if this leaf is in valid view space
[2017]2423                        BspLeaf *leaf = static_cast<BspLeaf *>(node);
[503]2424                        if (leaf->TreeValid() &&
[508]2425                                (!onlyUnmailed || !leaf->Mailed()) &&
[1707]2426                                ((maxPvsSize < 0) || (leaf->GetViewCell()->GetPvs().EvalPvsCost() <= maxPvsSize)))
[490]2427                        {
2428                                leaves.push_back(leaf);
2429                        }
[482]2430                }
2431                else
[463]2432                {
[2017]2433                        BspInterior *interior = static_cast<BspInterior *>(node);
[463]2434
2435                        nodeStack.push(interior->GetBack());
2436                        nodeStack.push(interior->GetFront());
2437                }
2438        }
2439}
2440
[489]2441
[463]2442AxisAlignedBox3 VspBspTree::GetBoundingBox() const
2443{
[1563]2444        return mBoundingBox;
[463]2445}
2446
[489]2447
[463]2448BspNode *VspBspTree::GetRoot() const
2449{
2450        return mRoot;
2451}
2452
[489]2453
[463]2454void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
2455{
2456        // the node became a leaf -> evaluate stats for leafs
[2017]2457        BspLeaf *leaf = static_cast<BspLeaf *>(data.mNode);
[463]2458
2459
[574]2460        if (data.mPvs > mBspStats.maxPvs)
[711]2461        {
[574]2462                mBspStats.maxPvs = data.mPvs;
[711]2463        }
2464
[656]2465        mBspStats.pvs += data.mPvs;
2466
[574]2467        if (data.mDepth < mBspStats.minDepth)
[711]2468        {
[574]2469                mBspStats.minDepth = data.mDepth;
[711]2470        }
[656]2471       
[463]2472        if (data.mDepth >= mTermMaxDepth)
[711]2473        {
[744]2474        ++ mBspStats.maxDepthNodes;
2475                //Debug << "new max depth: " << mBspStats.maxDepthNodes << endl;
[711]2476        }
[611]2477
[508]2478        // accumulate rays to compute rays /  leaf
[574]2479        mBspStats.accumRays += (int)data.mRays->size();
[463]2480
[437]2481        if (data.mPvs < mTermMinPvs)
[574]2482                ++ mBspStats.minPvsNodes;
[437]2483
2484        if ((int)data.mRays->size() < mTermMinRays)
[574]2485                ++ mBspStats.minRaysNodes;
[437]2486
2487        if (data.GetAvgRayContribution() > mTermMaxRayContribution)
[574]2488                ++ mBspStats.maxRayContribNodes;
[482]2489
[547]2490        if (data.mProbability <= mTermMinProbability)
[574]2491                ++ mBspStats.minProbabilityNodes;
[508]2492       
[474]2493        // accumulate depth to compute average depth
[574]2494        mBspStats.accumDepth += data.mDepth;
[463]2495
[612]2496        ++ mCreatedViewCells;
[656]2497
[1715]2498#ifdef GTP_DEBUG
[463]2499        Debug << "BSP stats: "
2500                  << "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
2501                  << "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
[1027]2502                  << "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
[1707]2503                  << "#pvs: " << leaf->GetViewCell()->GetPvs().EvalPvsCost() << "), "
[463]2504                  << "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
2505#endif
2506}
2507
[612]2508
[463]2509int VspBspTree::CastRay(Ray &ray)
2510{
2511        int hits = 0;
[482]2512
[600]2513        stack<BspRayTraversalData> tQueue;
[482]2514
[463]2515        float maxt, mint;
2516
[1563]2517        if (!mBoundingBox.GetRaySegment(ray, mint, maxt))
[463]2518                return 0;
2519
2520        Intersectable::NewMail();
[600]2521        ViewCell::NewMail();
[1012]2522
[463]2523        Vector3 entp = ray.Extrap(mint);
2524        Vector3 extp = ray.Extrap(maxt);
[482]2525
[463]2526        BspNode *node = mRoot;
2527        BspNode *farChild = NULL;
[482]2528
[463]2529        while (1)
2530        {
[482]2531                if (!node->IsLeaf())
[463]2532                {
[2017]2533                        BspInterior *in = static_cast<BspInterior *>(node);
[463]2534
2535                        Plane3 splitPlane = in->GetPlane();
2536                        const int entSide = splitPlane.Side(entp);
2537                        const int extSide = splitPlane.Side(extp);
2538
2539                        if (entSide < 0)
2540                        {
2541                                node = in->GetBack();
2542
2543                                if(extSide <= 0) // plane does not split ray => no far child
2544                                        continue;
[482]2545
[463]2546                                farChild = in->GetFront(); // plane splits ray
2547
2548                        } else if (entSide > 0)
2549                        {
2550                                node = in->GetFront();
2551
2552                                if (extSide >= 0) // plane does not split ray => no far child
2553                                        continue;
2554
[482]2555                                farChild = in->GetBack(); // plane splits ray
[463]2556                        }
2557                        else // ray and plane are coincident
2558                        {
[1012]2559                                // matt: WHAT TO DO IN THIS CASE ?
[463]2560                                //break;
2561                                node = in->GetFront();
2562                                continue;
2563                        }
2564
2565                        // push data for far child
[600]2566                        tQueue.push(BspRayTraversalData(farChild, extp, maxt));
[463]2567
2568                        // find intersection of ray segment with plane
2569                        float t;
2570                        extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
2571                        maxt *= t;
[1047]2572                }
2573                else // reached leaf => intersection with view cell
[463]2574                {
[2017]2575                        BspLeaf *leaf = static_cast<BspLeaf *>(node);
[482]2576
[463]2577                        if (!leaf->GetViewCell()->Mailed())
2578                        {
2579                                //ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
2580                                leaf->GetViewCell()->Mail();
2581                                ++ hits;
2582                        }
[482]2583
[463]2584                        //-- fetch the next far child from the stack
[600]2585                        if (tQueue.empty())
[463]2586                                break;
[482]2587
[463]2588                        entp = extp;
2589                        mint = maxt; // NOTE: need this?
2590
2591                        if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
2592                                break;
2593
[600]2594                        BspRayTraversalData &s = tQueue.top();
[463]2595
2596                        node = s.mNode;
2597                        extp = s.mExitPoint;
2598                        maxt = s.mMaxT;
2599
[600]2600                        tQueue.pop();
[463]2601                }
2602        }
2603
2604        return hits;
2605}
2606
[532]2607
[1072]2608void VspBspTree::CollectViewCells(ViewCellContainer &viewCells,
2609                                                                  bool onlyValid) const
[463]2610{
[532]2611        ViewCell::NewMail();
[551]2612        CollectViewCells(mRoot, onlyValid, viewCells, true);
[532]2613}
2614
2615
[1545]2616void VspBspTree::CollectViewCells(BspNode *root,
2617                                                                  bool onlyValid,
2618                                                                  ViewCellContainer &viewCells,
2619                                                                  bool onlyUnmailed) const
2620{
2621        stack<BspNode *> nodeStack;
2622
2623        if (!root)
2624                return;
2625
2626        nodeStack.push(root);
2627       
2628        while (!nodeStack.empty())
2629        {
2630                BspNode *node = nodeStack.top();
2631                nodeStack.pop();
2632               
2633                if (node->IsLeaf())
2634                {
2635                        if (!onlyValid || node->TreeValid())
2636                        {
[2017]2637                                ViewCellLeaf *leafVc = static_cast<BspLeaf *>(node)->GetViewCell();
[1545]2638
2639                                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leafVc);
2640                                               
2641                                if (!onlyUnmailed || !viewCell->Mailed())
2642                                {
2643                                        viewCell->Mail();
2644                                        viewCells.push_back(viewCell);
2645                                }
2646                        }
2647                }
2648                else
2649                {
[2017]2650                        BspInterior *interior = static_cast<BspInterior *>(node);
[1545]2651               
2652                        nodeStack.push(interior->GetFront());
2653                        nodeStack.push(interior->GetBack());
2654                }
2655        }
2656}
2657
2658
[574]2659void VspBspTree::CollapseViewCells()
[542]2660{
[590]2661// TODO
[728]2662#if HAS_TO_BE_REDONE
[542]2663        stack<BspNode *> nodeStack;
2664
2665        if (!mRoot)
2666                return;
2667
2668        nodeStack.push(mRoot);
2669       
2670        while (!nodeStack.empty())
2671        {
2672                BspNode *node = nodeStack.top();
2673                nodeStack.pop();
2674               
2675                if (node->IsLeaf())
[574]2676        {
[2017]2677                        BspViewCell *viewCell = static_cast<BspLeaf *>(node)->GetViewCell();
[542]2678
[574]2679                        if (!viewCell->GetValid())
[542]2680                        {
[2017]2681                                BspViewCell *viewCell = static_cast<BspLeaf *>(node)->GetViewCell();
[580]2682       
2683                                ViewCellContainer leaves;
[590]2684                                mViewCellsTree->CollectLeaves(viewCell, leaves);
[580]2685
2686                                ViewCellContainer::const_iterator it, it_end = leaves.end();
2687
2688                                for (it = leaves.begin(); it != it_end; ++ it)
[542]2689                                {
[2017]2690                                        BspLeaf *l = static_cast<BspViewCell *>(*it)->mLeaf;
[574]2691                                        l->SetViewCell(GetOrCreateOutOfBoundsCell());
2692                                        ++ mBspStats.invalidLeaves;
2693                                }
[542]2694
[574]2695                                // add to unbounded view cell
2696                                GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
2697                                DEL_PTR(viewCell);
2698                        }
2699                }
2700                else
2701                {
[2017]2702                        BspInterior *interior = static_cast<BspInterior *>(node);
[574]2703               
2704                        nodeStack.push(interior->GetFront());
2705                        nodeStack.push(interior->GetBack());
2706                }
2707        }
[542]2708
[574]2709        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
[590]2710#endif
[574]2711}
2712
2713
[639]2714void VspBspTree::CollectRays(VssRayContainer &rays)
2715{
2716        vector<BspLeaf *> leaves;
2717
2718        vector<BspLeaf *>::const_iterator lit, lit_end = leaves.end();
2719
2720        for (lit = leaves.begin(); lit != lit_end; ++ lit)
2721        {
2722                BspLeaf *leaf = *lit;
2723                VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
2724
2725                for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
2726                        rays.push_back(*rit);
2727        }
2728}
2729
2730
[574]2731void VspBspTree::ValidateTree()
2732{
2733        stack<BspNode *> nodeStack;
2734
2735        if (!mRoot)
2736                return;
2737
2738        nodeStack.push(mRoot);
2739       
2740        mBspStats.invalidLeaves = 0;
2741        while (!nodeStack.empty())
2742        {
2743                BspNode *node = nodeStack.top();
2744                nodeStack.pop();
2745               
2746                if (node->IsLeaf())
2747                {
[2017]2748                        BspLeaf *leaf = static_cast<BspLeaf *>(node);
[574]2749
2750                        if (!leaf->GetViewCell()->GetValid())
2751                                ++ mBspStats.invalidLeaves;
2752
2753                        // validity flags don't match => repair
2754                        if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
2755                        {
2756                                leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
2757                                PropagateUpValidity(leaf);
[542]2758                        }
2759                }
2760                else
2761                {
[2017]2762                        BspInterior *interior = static_cast<BspInterior *>(node);
[542]2763               
2764                        nodeStack.push(interior->GetFront());
2765                        nodeStack.push(interior->GetBack());
2766                }
2767        }
[562]2768
[574]2769        Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
[542]2770}
2771
[547]2772
[648]2773void VspBspTree::PreprocessPolygons(PolygonContainer &polys)
2774{
2775        // preprocess: throw out polygons coincident to the view space box (not needed)
2776        PolygonContainer boxPolys;
[1027]2777       
[1563]2778        mBoundingBox.ExtractPolys(boxPolys);
[648]2779        vector<Plane3> boxPlanes;
2780
2781        PolygonContainer::iterator pit, pit_end = boxPolys.end();
2782
2783        // extract planes of box
2784        // TODO: can be done more elegantly than first extracting polygons
2785        // and take their planes
2786        for (pit = boxPolys.begin(); pit != pit_end; ++ pit)
2787        {
2788                boxPlanes.push_back((*pit)->GetSupportingPlane());
2789        }
2790
2791        pit_end = polys.end();
2792
2793        for (pit = polys.begin(); pit != pit_end; ++ pit)
2794        {
2795                vector<Plane3>::const_iterator bit, bit_end = boxPlanes.end();
2796               
2797                for (bit = boxPlanes.begin(); (bit != bit_end) && (*pit); ++ bit)
2798                {
2799                        const int cf = (*pit)->ClassifyPlane(*bit, mEpsilon);
2800
2801                        if (cf == Polygon3::COINCIDENT)
2802                        {
2803                                DEL_PTR(*pit);
2804                                //Debug << "coincident!!" << endl;
2805                        }
2806                }
2807        }
2808
[1002]2809        // remove deleted entries after swapping them to end of vector
[648]2810        for (int i = 0; i < (int)polys.size(); ++ i)
2811        {
2812                while (!polys[i] && (i < (int)polys.size()))
2813                {
2814                        swap(polys[i], polys.back());
2815                        polys.pop_back();
2816                }
[1027]2817        }
[648]2818}
2819
2820
[463]2821float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
2822{
2823        float len = 0;
2824
2825        RayInfoContainer::const_iterator it, it_end = rays.end();
2826
2827        for (it = rays.begin(); it != it_end; ++ it)
2828                len += (*it).SegmentLength();
2829
2830        return len;
2831}
2832
[479]2833
[463]2834int VspBspTree::SplitRays(const Plane3 &plane,
[482]2835                                                  RayInfoContainer &rays,
2836                                                  RayInfoContainer &frontRays,
[639]2837                                                  RayInfoContainer &backRays) const
[463]2838{
2839        int splits = 0;
2840
[574]2841        RayInfoContainer::const_iterator it, it_end = rays.end();
2842
2843        for (it = rays.begin(); it != it_end; ++ it)
[463]2844        {
[574]2845                RayInfo bRay = *it;
2846               
[463]2847                VssRay *ray = bRay.mRay;
[473]2848                float t;
[463]2849
[485]2850                // get classification and receive new t
[463]2851                const int cf = bRay.ComputeRayIntersection(plane, t);
[482]2852
[463]2853                switch (cf)
2854                {
2855                case -1:
2856                        backRays.push_back(bRay);
2857                        break;
2858                case 1:
2859                        frontRays.push_back(bRay);
2860                        break;
[482]2861                case 0:
2862                        {
[485]2863                                //-- split ray
[639]2864                                //   test if start point behind or in front of plane
[485]2865                                const int side = plane.Side(bRay.ExtrapOrigin());
2866
[1011]2867                                ++ splits;
2868
[485]2869                                if (side <= 0)
2870                                {
2871                                        backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2872                                        frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2873                                }
2874                                else
2875                                {
2876                                        frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2877                                        backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2878                                }
[463]2879                        }
2880                        break;
2881                default:
[485]2882                        Debug << "Should not come here" << endl;
[463]2883                        break;
2884                }
2885        }
2886
2887        return splits;
2888}
2889
[479]2890
[463]2891void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
2892{
2893        BspNode *lastNode;
2894
2895        do
2896        {
2897                lastNode = n;
2898
2899                // want to get planes defining geometry of this node => don't take
2900                // split plane of node itself
2901                n = n->GetParent();
[482]2902
[463]2903                if (n)
2904                {
[2017]2905                        BspInterior *interior = static_cast<BspInterior *>(n);
2906                        Plane3 halfSpace = static_cast<BspInterior *>(interior)->GetPlane();
[463]2907
[683]2908            if (interior->GetBack() != lastNode)
[463]2909                                halfSpace.ReverseOrientation();
2910
2911                        halfSpaces.push_back(halfSpace);
2912                }
2913        }
2914        while (n);
2915}
2916
[485]2917
[482]2918void VspBspTree::ConstructGeometry(BspNode *n,
[503]2919                                                                   BspNodeGeometry &geom) const
[463]2920{
[437]2921        vector<Plane3> halfSpaces;
2922        ExtractHalfSpaces(n, halfSpaces);
2923
2924        PolygonContainer candidatePolys;
[678]2925        vector<Plane3> candidatePlanes;
[437]2926
[719]2927        vector<Plane3>::const_iterator pit, pit_end = halfSpaces.end();
2928
[683]2929        // bounded planes are added to the polygons
[719]2930        for (pit = halfSpaces.begin(); pit != pit_end; ++ pit)
[437]2931        {
[719]2932                Polygon3 *p = GetBoundingBox().CrossSection(*pit);
[482]2933
[448]2934                if (p->Valid(mEpsilon))
[437]2935                {
[683]2936                        candidatePolys.push_back(p);
[719]2937                        candidatePlanes.push_back(*pit);
[437]2938                }
2939        }
2940
2941        // add faces of bounding box (also could be faces of the cell)
2942        for (int i = 0; i < 6; ++ i)
2943        {
2944                VertexContainer vertices;
[482]2945
[437]2946                for (int j = 0; j < 4; ++ j)
[1563]2947                {
2948                        vertices.push_back(mBoundingBox.GetFace(i).mVertices[j]);
2949                }
[437]2950
[678]2951                Polygon3 *poly = new Polygon3(vertices);
2952
2953                candidatePolys.push_back(poly);
2954                candidatePlanes.push_back(poly->GetSupportingPlane());
[437]2955        }
2956
2957        for (int i = 0; i < (int)candidatePolys.size(); ++ i)
2958        {
2959                // polygon is split by all other planes
2960                for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
2961                {
2962                        if (i == j) // polygon and plane are coincident
2963                                continue;
2964
2965                        VertexContainer splitPts;
2966                        Polygon3 *frontPoly, *backPoly;
2967
[482]2968                        const int cf =
[448]2969                                candidatePolys[i]->ClassifyPlane(halfSpaces[j],
2970                                                                                                 mEpsilon);
[482]2971
[437]2972                        switch (cf)
2973                        {
2974                                case Polygon3::SPLIT:
2975                                        frontPoly = new Polygon3();
2976                                        backPoly = new Polygon3();
2977
[482]2978                                        candidatePolys[i]->Split(halfSpaces[j],
2979                                                                                         *frontPoly,
[448]2980                                                                                         *backPoly,
2981                                                                                         mEpsilon);
[437]2982
2983                                        DEL_PTR(candidatePolys[i]);
2984
[683]2985                                        if (backPoly->Valid(mEpsilon))
[1563]2986                                        {
[683]2987                                                candidatePolys[i] = backPoly;
[1563]2988                                        }
[437]2989                                        else
[1563]2990                                        {
[683]2991                                                DEL_PTR(backPoly);
[1563]2992                                        }
[437]2993
[683]2994                                        // outside, don't need this
2995                                        DEL_PTR(frontPoly);
[437]2996                                        break;
[1563]2997
[683]2998                                // polygon outside of halfspace
2999                                case Polygon3::FRONT_SIDE:
[437]3000                                        DEL_PTR(candidatePolys[i]);
3001                                        break;
[1563]3002
[437]3003                                // just take polygon as it is
[683]3004                                case Polygon3::BACK_SIDE:
[437]3005                                case Polygon3::COINCIDENT:
3006                                default:
3007                                        break;
3008                        }
3009                }
[482]3010
[437]3011                if (candidatePolys[i])
[678]3012                {
3013                        geom.Add(candidatePolys[i], candidatePlanes[i]);
3014                }
[437]3015        }
[463]3016}
3017
[485]3018
[1563]3019bool VspBspTree::IsOutOfBounds(ViewCell *vc) const
3020{
3021        return vc->GetId() == OUT_OF_BOUNDS_ID;
3022}
3023
3024
3025void VspBspTree::SetViewCellsTree(ViewCellsTree *viewCellsTree)
3026{
3027        mViewCellsTree = viewCellsTree;
3028}
3029
3030
[582]3031void VspBspTree::ConstructGeometry(ViewCell *vc,
[503]3032                                                                   BspNodeGeometry &vcGeom) const
[589]3033{
[1551]3034        // if false, cannot construct geometry for interior leaf
3035        if (!mViewCellsTree)
3036                return;
3037
[580]3038        ViewCellContainer leaves;
[590]3039        mViewCellsTree->CollectLeaves(vc, leaves);
[463]3040
[580]3041        ViewCellContainer::const_iterator it, it_end = leaves.end();
3042
[463]3043        for (it = leaves.begin(); it != it_end; ++ it)
[580]3044        {
[1563]3045                if (IsOutOfBounds(*it))
[1557]3046                        continue;
3047               
[2017]3048                BspViewCell *bspVc = static_cast<BspViewCell *>(*it);
[1551]3049                vector<BspLeaf *>::const_iterator bit, bit_end = bspVc->mLeaves.end();
3050
3051                for (bit = bspVc->mLeaves.begin(); bit != bit_end; ++ bit)
3052                {
3053                        BspLeaf *l = *bit;
3054                        ConstructGeometry(l, vcGeom);
3055                }
[580]3056        }
[463]3057}
3058
[485]3059
[482]3060int VspBspTree::FindNeighbors(BspNode *n, vector<BspLeaf *> &neighbors,
[562]3061                                                          const bool onlyUnmailed) const
[463]3062{
[551]3063        stack<bspNodePair> nodeStack;
3064       
3065        BspNodeGeometry nodeGeom;
3066        ConstructGeometry(n, nodeGeom);
[801]3067//      const float eps = 0.5f;
[752]3068        const float eps = 0.01f;
[500]3069        // split planes from the root to this node
3070        // needed to verify that we found neighbor leaf
[557]3071        // TODO: really needed?
[463]3072        vector<Plane3> halfSpaces;
3073        ExtractHalfSpaces(n, halfSpaces);
3074
[551]3075
3076        BspNodeGeometry *rgeom = new BspNodeGeometry();
3077        ConstructGeometry(mRoot, *rgeom);
3078
3079        nodeStack.push(bspNodePair(mRoot, rgeom));
3080
[482]3081        while (!nodeStack.empty())
[463]3082        {
[551]3083                BspNode *node = nodeStack.top().first;
3084                BspNodeGeometry *geom = nodeStack.top().second;
[562]3085       
[463]3086                nodeStack.pop();
3087
[557]3088                if (node->IsLeaf())
[562]3089                {
[557]3090                        // test if this leaf is in valid view space
3091                        if (node->TreeValid() &&
3092                                (node != n) &&
3093                                (!onlyUnmailed || !node->Mailed()))
3094                        {
3095                                bool isAdjacent = true;
[551]3096
[570]3097                                if (1)
[557]3098                                {
[562]3099                                        // test all planes of current node if still adjacent
3100                                        for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
3101                                        {
3102                                                const int cf =
[678]3103                                                        Polygon3::ClassifyPlane(geom->GetPolys(),
[562]3104                                                                                                        halfSpaces[i],
[752]3105                                                                                                        eps);
[482]3106
[683]3107                                                if (cf == Polygon3::FRONT_SIDE)
[562]3108                                                {
3109                                                        isAdjacent = false;
3110                                                }
[557]3111                                        }
3112                                }
[562]3113                                else
[557]3114                                {
[562]3115                                        // TODO: why is this wrong??
3116                                        // test all planes of current node if still adjacent
[678]3117                                        for (int i = 0; (i < nodeGeom.Size()) && isAdjacent; ++ i)
[562]3118                                        {
[678]3119                                                Polygon3 *poly = nodeGeom.GetPolys()[i];
[555]3120
[562]3121                                                const int cf =
[678]3122                                                        Polygon3::ClassifyPlane(geom->GetPolys(),
[562]3123                                                                                                        poly->GetSupportingPlane(),
[752]3124                                                                                                        eps);
[557]3125
[683]3126                                                if (cf == Polygon3::FRONT_SIDE)
[562]3127                                                {
3128                                                        isAdjacent = false;
3129                                                }
[557]3130                                        }
[570]3131                                }
[557]3132                                // neighbor was found
3133                                if (isAdjacent)
[562]3134                                {       
[2017]3135                                        neighbors.push_back(static_cast<BspLeaf *>(node));
[562]3136                                }
[463]3137                        }
[562]3138                }
3139                else
3140                {
[2017]3141                        BspInterior *interior = static_cast<BspInterior *>(node);
[482]3142
[678]3143                        const int cf = Polygon3::ClassifyPlane(nodeGeom.GetPolys(),
[562]3144                                                                                                   interior->GetPlane(),
[752]3145                                                                                                   eps);
[551]3146                       
[562]3147                        BspNode *front = interior->GetFront();
3148                        BspNode *back = interior->GetBack();
[551]3149           
[562]3150                        BspNodeGeometry *fGeom = new BspNodeGeometry();
3151                        BspNodeGeometry *bGeom = new BspNodeGeometry();
[463]3152
[562]3153                        geom->SplitGeometry(*fGeom,
3154                                                                *bGeom,
3155                                                                interior->GetPlane(),
[1563]3156                                                                mBoundingBox,
[675]3157                                                                //0.0000001f);
[752]3158                                                                eps);
[551]3159               
[683]3160                        if (cf == Polygon3::BACK_SIDE)
[562]3161                        {
[683]3162                                nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
3163                                DEL_PTR(fGeom);
[562]3164                        }
3165                        else
3166                        {
[683]3167                                if (cf == Polygon3::FRONT_SIDE)
[551]3168                                {
[683]3169                                        nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
3170                                        DEL_PTR(bGeom);
[551]3171                                }
[482]3172                                else
[562]3173                                {       // random decision
3174                                        nodeStack.push(bspNodePair(front, fGeom));
3175                                        nodeStack.push(bspNodePair(back, bGeom));
[463]3176                                }
[551]3177                        }
[463]3178                }
[562]3179       
[551]3180                DEL_PTR(geom);
[463]3181        }
[482]3182
[463]3183        return (int)neighbors.size();
3184}
3185
[489]3186
[600]3187
[710]3188int VspBspTree::FindApproximateNeighbors(BspNode *n,
3189                                                                                 vector<BspLeaf *> &neighbors,
[600]3190                                                                                 const bool onlyUnmailed) const
3191{
3192        stack<bspNodePair> nodeStack;
3193       
3194        BspNodeGeometry nodeGeom;
3195        ConstructGeometry(n, nodeGeom);
3196       
[752]3197        float eps = 0.01f;
[600]3198        // split planes from the root to this node
3199        // needed to verify that we found neighbor leaf
3200        // TODO: really needed?
3201        vector<Plane3> halfSpaces;
3202        ExtractHalfSpaces(n, halfSpaces);
3203
3204
3205        BspNodeGeometry *rgeom = new BspNodeGeometry();
3206        ConstructGeometry(mRoot, *rgeom);
3207
3208        nodeStack.push(bspNodePair(mRoot, rgeom));
3209
3210        while (!nodeStack.empty())
3211        {
3212                BspNode *node = nodeStack.top().first;
3213                BspNodeGeometry *geom = nodeStack.top().second;
3214       
3215                nodeStack.pop();
3216
3217                if (node->IsLeaf())
3218                {
3219                        // test if this leaf is in valid view space
3220                        if (node->TreeValid() &&
3221                                (node != n) &&
3222                                (!onlyUnmailed || !node->Mailed()))
3223                        {
3224                                bool isAdjacent = true;
3225
3226                                // neighbor was found
3227                                if (isAdjacent)
3228                                {       
[2017]3229                                        neighbors.push_back(static_cast<BspLeaf *>(node));
[600]3230                                }
3231                        }
3232                }
3233                else
3234                {
[2017]3235                        BspInterior *interior = static_cast<BspInterior *>(node);
[600]3236
[678]3237                        const int cf = Polygon3::ClassifyPlane(nodeGeom.GetPolys(),
[600]3238                                                                                                   interior->GetPlane(),
[752]3239                                                                                                   eps);
[600]3240                       
3241                        BspNode *front = interior->GetFront();
3242                        BspNode *back = interior->GetBack();
3243           
3244                        BspNodeGeometry *fGeom = new BspNodeGeometry();
3245                        BspNodeGeometry *bGeom = new BspNodeGeometry();
3246
3247                        geom->SplitGeometry(*fGeom,
3248                                                                *bGeom,
3249                                                                interior->GetPlane(),
[1563]3250                                                                mBoundingBox,
[675]3251                                                                //0.0000001f);
[752]3252                                                                eps);
[600]3253               
[683]3254                        if (cf == Polygon3::BACK_SIDE)
[600]3255                        {
[683]3256                                nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
3257                                DEL_PTR(fGeom);
[710]3258                                }
[600]3259                        else
3260                        {
[683]3261                                if (cf == Polygon3::FRONT_SIDE)
[600]3262                                {
[683]3263                                        nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
3264                                        DEL_PTR(bGeom);
[600]3265                                }
3266                                else
3267                                {       // random decision
3268                                        nodeStack.push(bspNodePair(front, fGeom));
3269                                        nodeStack.push(bspNodePair(back, bGeom));
3270                                }
3271                        }
3272                }
3273       
3274                DEL_PTR(geom);
3275        }
3276
3277        return (int)neighbors.size();
3278}
3279
3280
3281
[463]3282BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
3283{
3284    stack<BspNode *> nodeStack;
3285        nodeStack.push(mRoot);
[482]3286
[463]3287        int mask = rand();
[482]3288
3289        while (!nodeStack.empty())
[463]3290        {
3291                BspNode *node = nodeStack.top();
3292                nodeStack.pop();
[482]3293
3294                if (node->IsLeaf())
[463]3295                {
[2017]3296                        return static_cast<BspLeaf *>(node);
[482]3297                }
3298                else
[463]3299                {
[2017]3300                        BspInterior *interior = static_cast<BspInterior *>(node);
[463]3301                        BspNode *next;
[503]3302                        BspNodeGeometry geom;
[482]3303
[463]3304                        // todo: not very efficient: constructs full cell everytime
[498]3305                        ConstructGeometry(interior, geom);
[463]3306
[503]3307                        const int cf =
[678]3308                                Polygon3::ClassifyPlane(geom.GetPolys(), halfspace, mEpsilon);
[463]3309
3310                        if (cf == Polygon3::BACK_SIDE)
3311                                next = interior->GetFront();
3312                        else
3313                                if (cf == Polygon3::FRONT_SIDE)
3314                                        next = interior->GetFront();
[482]3315                        else
[463]3316                        {
3317                                // random decision
3318                                if (mask & 1)
3319                                        next = interior->GetBack();
3320                                else
3321                                        next = interior->GetFront();
3322                                mask = mask >> 1;
3323                        }
3324
3325                        nodeStack.push(next);
3326                }
3327        }
[482]3328
[463]3329        return NULL;
3330}
3331
[694]3332
[463]3333BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
3334{
3335        stack<BspNode *> nodeStack;
[482]3336
[463]3337        nodeStack.push(mRoot);
3338
3339        int mask = rand();
[482]3340
3341        while (!nodeStack.empty())
[463]3342        {
3343                BspNode *node = nodeStack.top();
3344                nodeStack.pop();
[482]3345
3346                if (node->IsLeaf())
[463]3347                {
3348                        if ( (!onlyUnmailed || !node->Mailed()) )
[2017]3349                                return static_cast<BspLeaf *>(node);
[463]3350                }
[482]3351                else
[463]3352                {
[2017]3353                        BspInterior *interior = static_cast<BspInterior *>(node);
[463]3354
3355                        // random decision
3356                        if (mask & 1)
3357                                nodeStack.push(interior->GetBack());
3358                        else
3359                                nodeStack.push(interior->GetFront());
3360
3361                        mask = mask >> 1;
3362                }
3363        }
[482]3364
[463]3365        return NULL;
3366}
3367
[694]3368
[463]3369int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
3370{
3371        int pvsSize = 0;
3372
3373        RayInfoContainer::const_iterator rit, rit_end = rays.end();
3374
3375        Intersectable::NewMail();
3376
3377        for (rit = rays.begin(); rit != rays.end(); ++ rit)
3378        {
3379                VssRay *ray = (*rit).mRay;
[482]3380
[1020]3381                if (COUNT_ORIGIN_OBJECTS && ray->mOriginObject)
[463]3382                {
3383                        if (!ray->mOriginObject->Mailed())
3384                        {
3385                                ray->mOriginObject->Mail();
3386                                ++ pvsSize;
3387                        }
3388                }
[1020]3389
[463]3390                if (ray->mTerminationObject)
3391                {
3392                        if (!ray->mTerminationObject->Mailed())
3393                        {
3394                                ray->mTerminationObject->Mail();
3395                                ++ pvsSize;
3396                        }
3397                }
3398        }
3399
3400        return pvsSize;
3401}
3402
[694]3403
[463]3404float VspBspTree::GetEpsilon() const
3405{
3406        return mEpsilon;
3407}
3408
3409
3410int VspBspTree::SplitPolygons(const Plane3 &plane,
[482]3411                                                          PolygonContainer &polys,
3412                                                          PolygonContainer &frontPolys,
3413                                                          PolygonContainer &backPolys,
[463]3414                                                          PolygonContainer &coincident) const
3415{
3416        int splits = 0;
3417
[574]3418        PolygonContainer::const_iterator it, it_end = polys.end();
3419
3420        for (it = polys.begin(); it != polys.end(); ++ it)     
[463]3421        {
[574]3422                Polygon3 *poly = *it;
[463]3423
3424                // classify polygon
3425                const int cf = poly->ClassifyPlane(plane, mEpsilon);
3426
3427                switch (cf)
3428                {
3429                        case Polygon3::COINCIDENT:
3430                                coincident.push_back(poly);
[482]3431                                break;
3432                        case Polygon3::FRONT_SIDE:
[463]3433                                frontPolys.push_back(poly);
3434                                break;
3435                        case Polygon3::BACK_SIDE:
3436                                backPolys.push_back(poly);
3437                                break;
3438                        case Polygon3::SPLIT:
3439                                backPolys.push_back(poly);
3440                                frontPolys.push_back(poly);
3441                                ++ splits;
3442                                break;
3443                        default:
3444                Debug << "SHOULD NEVER COME HERE\n";
3445                                break;
3446                }
3447        }
3448
3449        return splits;
3450}
[466]3451
3452
[469]3453int VspBspTree::CastLineSegment(const Vector3 &origin,
3454                                                                const Vector3 &termination,
[882]3455                                                                ViewCellContainer &viewcells)
[466]3456{
[469]3457        int hits = 0;
[719]3458        stack<BspRayTraversalData> tStack;
[482]3459
[469]3460        float mint = 0.0f, maxt = 1.0f;
[482]3461
[1291]3462        //ViewCell::NewMail();
[482]3463
[469]3464        Vector3 entp = origin;
3465        Vector3 extp = termination;
[482]3466
[469]3467        BspNode *node = mRoot;
3468        BspNode *farChild = NULL;
[482]3469
[485]3470        float t;
[752]3471        const float thresh = 1e-6f; // matt: change this to adjustable value
[2015]3472
[482]3473        while (1)
[469]3474        {
[482]3475                if (!node->IsLeaf())
[469]3476                {
[2017]3477                        BspInterior *in = static_cast<BspInterior *>(node);
[482]3478
[2021]3479                        Plane3 &splitPlane = in->GetPlane();
[485]3480                       
[666]3481                        const int entSide = splitPlane.Side(entp, thresh);
3482                        const int extSide = splitPlane.Side(extp, thresh);
[482]3483
[485]3484                        if (entSide < 0)
[469]3485                        {
[666]3486                                node = in->GetBack();
3487                               
3488                                // plane does not split ray => no far child
3489                                if (extSide <= 0)
3490                                        continue;
3491 
3492                                farChild = in->GetFront(); // plane splits ray
[485]3493                        }
3494                        else if (entSide > 0)
[469]3495                        {
[666]3496                                node = in->GetFront();
[482]3497
[666]3498                                if (extSide >= 0) // plane does not split ray => no far child
3499                                        continue;
[482]3500
3501                                farChild = in->GetBack(); // plane splits ray
[469]3502                        }
[694]3503                        else // one of the ray end points is on the plane
[1020]3504                        {       
3505                                // NOTE: what to do if ray is coincident with plane?
[485]3506                                if (extSide < 0)
3507                                        node = in->GetBack();
[694]3508                                else //if (extSide > 0)
[485]3509                                        node = in->GetFront();
[694]3510                                //else break; // coincident => count no intersections
3511
[485]3512                                continue; // no far child
[469]3513                        }
[482]3514
[469]3515                        // push data for far child
[719]3516                        tStack.push(BspRayTraversalData(farChild, extp));
[482]3517
[469]3518                        // find intersection of ray segment with plane
3519                        extp = splitPlane.FindIntersection(origin, extp, &t);
[485]3520                }
3521                else
[469]3522                {
3523                        // reached leaf => intersection with view cell
[2017]3524                        BspLeaf *leaf = static_cast<BspLeaf *>(node);
[666]3525                        ViewCell *viewCell;
3526                       
[1020]3527                        // question: always contribute to leaf or to currently active view cell?
[2021]3528                        //                      if (0)
3529                        //                              viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
3530                        //                      else
3531                        viewCell = leaf->GetViewCell();
[482]3532
[590]3533                        if (!viewCell->Mailed())
[469]3534                        {
[590]3535                                viewcells.push_back(viewCell);
3536                                viewCell->Mail();
[469]3537                                ++ hits;
3538                        }
[482]3539
[469]3540                        //-- fetch the next far child from the stack
[719]3541                        if (tStack.empty())
[469]3542                                break;
[482]3543
[469]3544                        entp = extp;
[485]3545                       
[719]3546                        const BspRayTraversalData &s = tStack.top();
[482]3547
[469]3548                        node = s.mNode;
3549                        extp = s.mExitPoint;
[482]3550
[719]3551                        tStack.pop();
[469]3552                }
[466]3553        }
[487]3554
[469]3555        return hits;
[466]3556}
[478]3557
[576]3558
3559
3560
[485]3561int VspBspTree::TreeDistance(BspNode *n1, BspNode *n2) const
[482]3562{
3563        std::deque<BspNode *> path1;
3564        BspNode *p1 = n1;
[479]3565
[482]3566        // create path from node 1 to root
3567        while (p1)
3568        {
3569                if (p1 == n2) // second node on path
3570                        return (int)path1.size();
3571
3572                path1.push_front(p1);
3573                p1 = p1->GetParent();
3574        }
3575
3576        int depth = n2->GetDepth();
3577        int d = depth;
3578
3579        BspNode *p2 = n2;
3580
3581        // compare with same depth
3582        while (1)
3583        {
3584                if ((d < (int)path1.size()) && (p2 == path1[d]))
3585                        return (depth - d) + ((int)path1.size() - 1 - d);
3586
3587                -- d;
3588                p2 = p2->GetParent();
3589        }
3590
3591        return 0; // never come here
3592}
3593
[580]3594
[501]3595BspNode *VspBspTree::CollapseTree(BspNode *node, int &collapsed)
[479]3596{
[590]3597// TODO
[728]3598#if HAS_TO_BE_REDONE
[495]3599        if (node->IsLeaf())
[479]3600                return node;
3601
[2017]3602        BspInterior *interior = static_cast<BspInterior *>(node);
[492]3603
[501]3604        BspNode *front = CollapseTree(interior->GetFront(), collapsed);
3605        BspNode *back = CollapseTree(interior->GetBack(), collapsed);
[492]3606
[479]3607        if (front->IsLeaf() && back->IsLeaf())
3608        {
[2017]3609                BspLeaf *frontLeaf = static_cast<BspLeaf *>(front);
3610                BspLeaf *backLeaf = static_cast<BspLeaf *>(back);
[479]3611
3612                //-- collapse tree
3613                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
3614                {
3615                        BspViewCell *vc = frontLeaf->GetViewCell();
3616
3617                        BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);
[489]3618                        leaf->SetTreeValid(frontLeaf->TreeValid());
[482]3619
[479]3620                        // replace a link from node's parent
3621                        if (leaf->GetParent())
3622                                leaf->GetParent()->ReplaceChildLink(node, leaf);
[517]3623                        else
3624                                mRoot = leaf;
3625
[501]3626                        ++ collapsed;
[479]3627                        delete interior;
3628
3629                        return leaf;
3630                }
3631        }
[590]3632#endif
[495]3633        return node;
3634}
3635
3636
[501]3637int VspBspTree::CollapseTree()
[495]3638{
[501]3639        int collapsed = 0;
[580]3640        //TODO
[728]3641#if HAS_TO_BE_REDONE
[501]3642        (void)CollapseTree(mRoot, collapsed);
[517]3643
[485]3644        // revalidate leaves
[517]3645        RepairViewCellsLeafLists();
[580]3646#endif
[501]3647        return collapsed;
[479]3648}
3649
3650
[517]3651void VspBspTree::RepairViewCellsLeafLists()
[492]3652{
[590]3653// TODO
[728]3654#if HAS_TO_BE_REDONE
[479]3655        // list not valid anymore => clear
[492]3656        stack<BspNode *> nodeStack;
3657        nodeStack.push(mRoot);
3658
3659        ViewCell::NewMail();
3660
3661        while (!nodeStack.empty())
3662        {
3663                BspNode *node = nodeStack.top();
3664                nodeStack.pop();
3665
3666                if (node->IsLeaf())
3667                {
[2017]3668                        BspLeaf *leaf = static_cast<BspLeaf *>(node);
[492]3669
3670                        BspViewCell *viewCell = leaf->GetViewCell();
[590]3671
[492]3672                        if (!viewCell->Mailed())
3673                        {
3674                                viewCell->mLeaves.clear();
3675                                viewCell->Mail();
3676                        }
[580]3677       
[492]3678                        viewCell->mLeaves.push_back(leaf);
[590]3679
[492]3680                }
3681                else
3682                {
[2017]3683                        BspInterior *interior = static_cast<BspInterior *>(node);
[492]3684
3685                        nodeStack.push(interior->GetFront());
3686                        nodeStack.push(interior->GetBack());
3687                }
[479]3688        }
[590]3689// TODO
3690#endif
[479]3691}
3692
3693
[532]3694int VspBspTree::CastBeam(Beam &beam)
3695{
3696    stack<bspNodePair> nodeStack;
3697        BspNodeGeometry *rgeom = new BspNodeGeometry();
3698        ConstructGeometry(mRoot, *rgeom);
3699
3700        nodeStack.push(bspNodePair(mRoot, rgeom));
3701 
3702        ViewCell::NewMail();
3703
3704        while (!nodeStack.empty())
3705        {
3706                BspNode *node = nodeStack.top().first;
3707                BspNodeGeometry *geom = nodeStack.top().second;
3708                nodeStack.pop();
3709               
3710                AxisAlignedBox3 box;
[697]3711                geom->GetBoundingBox(box);
[532]3712
[535]3713                const int side = beam.ComputeIntersection(box);
[532]3714               
3715                switch (side)
3716                {
3717                case -1:
[547]3718                        CollectViewCells(node, true, beam.mViewCells, true);
[532]3719                        break;
3720                case 0:
[535]3721                       
[532]3722                        if (node->IsLeaf())
3723                        {
[2017]3724                                BspLeaf *leaf = static_cast<BspLeaf *>(node);
[535]3725                       
[532]3726                                if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
[535]3727                                {
3728                                        leaf->GetViewCell()->Mail();
[532]3729                                        beam.mViewCells.push_back(leaf->GetViewCell());
[535]3730                                }
[532]3731                        }
3732                        else
3733                        {
[2017]3734                                BspInterior *interior = static_cast<BspInterior *>(node);
[535]3735                       
[538]3736                                BspNode *first = interior->GetFront();
3737                                BspNode *second = interior->GetBack();
[535]3738           
3739                                BspNodeGeometry *firstGeom = new BspNodeGeometry();
3740                                BspNodeGeometry *secondGeom = new BspNodeGeometry();
3741
[538]3742                                geom->SplitGeometry(*firstGeom,
3743                                                                        *secondGeom,
3744                                                                        interior->GetPlane(),
[1563]3745                                                                        mBoundingBox,
[675]3746                                                                        //0.0000001f);
[538]3747                                                                        mEpsilon);
[535]3748
[532]3749                                // decide on the order of the nodes
3750                                if (DotProd(beam.mPlanes[0].mNormal,
3751                                        interior->GetPlane().mNormal) > 0)
3752                                {
3753                                        swap(first, second);
[535]3754                                        swap(firstGeom, secondGeom);
[532]3755                                }
3756
[535]3757                                nodeStack.push(bspNodePair(first, firstGeom));
3758                                nodeStack.push(bspNodePair(second, secondGeom));
[532]3759                        }
[535]3760                       
[532]3761                        break;
[538]3762                default:
[532]3763                        // default: cull
[538]3764                        break;
[532]3765                }
[538]3766               
[532]3767                DEL_PTR(geom);
[535]3768               
[532]3769        }
3770
[538]3771        return (int)beam.mViewCells.size();
[532]3772}
3773
3774
[485]3775void VspBspTree::SetViewCellsManager(ViewCellsManager *vcm)
[478]3776{
[485]3777        mViewCellsManager = vcm;
3778}
3779
3780
[580]3781int VspBspTree::CollectMergeCandidates(const vector<BspLeaf *> leaves,
3782                                                                           vector<MergeCandidate> &candidates)
[485]3783{
[478]3784        BspLeaf::NewMail();
[508]3785       
[478]3786        vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
3787
[580]3788        int numCandidates = 0;
[509]3789
[478]3790        // find merge candidates and push them into queue
3791        for (it = leaves.begin(); it != it_end; ++ it)
3792        {
3793                BspLeaf *leaf = *it;
[485]3794               
3795                // the same leaves must not be part of two merge candidates
3796                leaf->Mail();
[650]3797               
[485]3798                vector<BspLeaf *> neighbors;
[710]3799               
3800                // appoximate neighbor search has slightl relaxed constraints
[694]3801                if (1)
[650]3802                        FindNeighbors(leaf, neighbors, true);
3803                else
3804                        FindApproximateNeighbors(leaf, neighbors, true);
[710]3805
[485]3806                vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();
3807
3808                // TODO: test if at least one ray goes from one leaf to the other
3809                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
[508]3810                {
3811                        if ((*nit)->GetViewCell() != leaf->GetViewCell())
[509]3812                        {
[580]3813                                MergeCandidate mc(leaf->GetViewCell(), (*nit)->GetViewCell());
[564]3814
[1006]3815                                if (!leaf->GetViewCell()->GetPvs().Empty() ||
[710]3816                                        !(*nit)->GetViewCell()->GetPvs().Empty() ||
[676]3817                    leaf->IsSibling(*nit))
3818                                {
3819                                        candidates.push_back(mc);
3820                                }
3821
[580]3822                                ++ numCandidates;
3823                                if ((numCandidates % 1000) == 0)
[566]3824                                {
[580]3825                                        cout << "collected " << numCandidates << " merge candidates" << endl;
[566]3826                                }
[509]3827                        }
[485]3828                }
3829        }
3830
[580]3831        Debug << "merge queue: " << (int)candidates.size() << endl;
3832        Debug << "leaves in queue: " << numCandidates << endl;
3833       
[508]3834
[485]3835        return (int)leaves.size();
3836}
3837
3838
[580]3839int VspBspTree::CollectMergeCandidates(const VssRayContainer &rays,
3840                                                                           vector<MergeCandidate> &candidates)
[485]3841{
[547]3842        ViewCell::NewMail();
[503]3843        long startTime = GetTime();
[574]3844       
[485]3845        map<BspLeaf *, vector<BspLeaf*> > neighborMap;
[574]3846        ViewCellContainer::const_iterator iit;
[485]3847
[503]3848        int numLeaves = 0;
[485]3849       
3850        BspLeaf::NewMail();
3851
[574]3852        for (int i = 0; i < (int)rays.size(); ++ i)
[485]3853        { 
[574]3854                VssRay *ray = rays[i];
[547]3855       
[485]3856                // traverse leaves stored in the rays and compare and
3857                // merge consecutive leaves (i.e., the neighbors in the tree)
[574]3858                if (ray->mViewCells.size() < 2)
[485]3859                        continue;
[1551]3860
[574]3861                iit = ray->mViewCells.begin();
[2017]3862                BspViewCell *bspVc = static_cast<BspViewCell *>(*(iit ++));
[1551]3863                BspLeaf *leaf = bspVc->mLeaves[0];
[485]3864               
3865                // traverse intersections
[489]3866                // consecutive leaves are neighbors => add them to queue
[574]3867                for (; iit != ray->mViewCells.end(); ++ iit)
[485]3868                {
[489]3869                        // next pair
3870                        BspLeaf *prevLeaf = leaf;
[2017]3871                        bspVc = static_cast<BspViewCell *>(*iit);
[1551]3872            leaf = bspVc->mLeaves[0]; // exactly one leaf
[489]3873
[508]3874                        // view space not valid or same view cell
3875                        if (!leaf->TreeValid() || !prevLeaf->TreeValid() ||
3876                                (leaf->GetViewCell() == prevLeaf->GetViewCell()))
[489]3877                                continue;
3878
[580]3879                vector<BspLeaf *> &neighbors = neighborMap[leaf];
[485]3880                       
3881                        bool found = false;
[478]3882
[485]3883                        // both leaves inserted in queue already =>
3884                        // look if double pair already exists
3885                        if (leaf->Mailed() && prevLeaf->Mailed())
[478]3886                        {
[485]3887                                vector<BspLeaf *>::const_iterator it, it_end = neighbors.end();
3888                               
3889                for (it = neighbors.begin(); !found && (it != it_end); ++ it)
3890                                        if (*it == prevLeaf)
3891                                                found = true; // already in queue
3892                        }
[547]3893               
[485]3894                        if (!found)
3895                        {
[564]3896                                // this pair is not in map yet
[485]3897                                // => insert into the neighbor map and the queue
3898                                neighbors.push_back(prevLeaf);
3899                                neighborMap[prevLeaf].push_back(leaf);
[478]3900
[485]3901                                leaf->Mail();
3902                                prevLeaf->Mail();
[547]3903               
[580]3904                                MergeCandidate mc(leaf->GetViewCell(), prevLeaf->GetViewCell());
3905                               
3906                                candidates.push_back(mc);
[564]3907
[580]3908                                if (((int)candidates.size() % 1000) == 0)
[564]3909                                {
[580]3910                                        cout << "collected " << (int)candidates.size() << " merge candidates" << endl;
[564]3911                                }
[478]3912                        }
[485]3913        }
3914        }
[564]3915
[485]3916        Debug << "neighbormap size: " << (int)neighborMap.size() << endl;
[580]3917        Debug << "merge queue: " << (int)candidates.size() << endl;
[503]3918        Debug << "leaves in queue: " << numLeaves << endl;
[485]3919
[580]3920
[503]3921        //-- collect the leaves which haven't been found by ray casting
[542]3922        if (0)
3923        {
[551]3924                cout << "finding additional merge candidates using geometry" << endl;
[542]3925                vector<BspLeaf *> leaves;
[547]3926                CollectLeaves(leaves, true);
[542]3927                Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
[580]3928                CollectMergeCandidates(leaves, candidates);
[542]3929        }
[503]3930
3931        return numLeaves;
[485]3932}
3933
3934
3935
3936
[879]3937ViewCell *VspBspTree::GetViewCell(const Vector3 &point, const bool active)
[492]3938{
[879]3939        if (mRoot == NULL)
3940                return NULL;
3941
3942        stack<BspNode *> nodeStack;
3943        nodeStack.push(mRoot);
[492]3944 
[882]3945        ViewCellLeaf *viewcell = NULL;
[492]3946 
[879]3947        while (!nodeStack.empty()) 
3948        {
3949                BspNode *node = nodeStack.top();
3950                nodeStack.pop();
[492]3951       
[879]3952                if (node->IsLeaf())
3953                {
[2017]3954                        viewcell = static_cast<BspLeaf *>(node)->GetViewCell();
[879]3955                        break;
3956                }
3957                else   
3958                {       
[2017]3959                        BspInterior *interior = static_cast<BspInterior *>(node);
[879]3960       
3961                        // random decision
3962                        if (interior->GetPlane().Side(point) < 0)
3963                                nodeStack.push(interior->GetBack());
3964                        else
3965                                nodeStack.push(interior->GetFront());
3966                }
[492]3967        }
3968 
[879]3969        if (active)
3970                return mViewCellsTree->GetActiveViewCell(viewcell);
3971        else
3972                return viewcell;
[492]3973}
3974
3975
[487]3976bool VspBspTree::ViewPointValid(const Vector3 &viewPoint) const
3977{
3978        BspNode *node = mRoot;
[485]3979
[487]3980        while (1)
3981        {
3982                // early exit
3983                if (node->TreeValid())
3984                        return true;
3985
3986                if (node->IsLeaf())
3987                        return false;
3988                       
[2017]3989                BspInterior *in = static_cast<BspInterior *>(node);
[490]3990                                       
3991                if (in->GetPlane().Side(viewPoint) <= 0)
[487]3992                {
3993                        node = in->GetBack();
3994                }
3995                else
3996                {
3997                        node = in->GetFront();
3998                }
3999        }
4000
4001        // should never come here
4002        return false;
4003}
4004
4005
4006void VspBspTree::PropagateUpValidity(BspNode *node)
4007{
[574]4008        const bool isValid = node->TreeValid();
4009
4010        // propagative up invalid flag until only invalid nodes exist over this node
4011        if (!isValid)
[487]4012        {
[574]4013                while (!node->IsRoot() && node->GetParent()->TreeValid())
4014                {
4015                        node = node->GetParent();
4016                        node->SetTreeValid(false);
4017                }
[487]4018        }
[574]4019        else
4020        {
4021                // propagative up valid flag until one of the subtrees is invalid
4022                while (!node->IsRoot() && !node->TreeValid())
4023                {
4024            node = node->GetParent();
[2017]4025                        BspInterior *interior = static_cast<BspInterior *>(node);
[574]4026                       
4027                        // the parent is valid iff both leaves are valid
4028                        node->SetTreeValid(interior->GetBack()->TreeValid() &&
4029                                                           interior->GetFront()->TreeValid());
4030                }
4031        }
[487]4032}
4033
[1201]4034
4035bool VspBspTree::Export(OUT_STREAM &stream)
[503]4036{
[508]4037        ExportNode(mRoot, stream);
[503]4038        return true;
4039}
4040
[1201]4041
4042void VspBspTree::ExportNode(BspNode *node, OUT_STREAM &stream)
[508]4043{
4044        if (node->IsLeaf())
[503]4045        {
[2017]4046                BspLeaf *leaf = static_cast<BspLeaf *>(node);
[590]4047                ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
4048
[508]4049                int id = -1;
[590]4050                if (viewCell != mOutOfBoundsCell)
4051                        id = viewCell->GetId();
[503]4052
[508]4053                stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
[503]4054        }
[508]4055        else
[503]4056        {
[2017]4057                BspInterior *interior = static_cast<BspInterior *>(node);
[508]4058       
4059                Plane3 plane = interior->GetPlane();
4060                stream << "<Interior plane=\"" << plane.mNormal.x << " "
4061                           << plane.mNormal.y << " " << plane.mNormal.z << " "
4062                           << plane.mD << "\">" << endl;
[503]4063
[508]4064                ExportNode(interior->GetBack(), stream);
4065                ExportNode(interior->GetFront(), stream);
[503]4066
[508]4067                stream << "</Interior>" << endl;
[503]4068        }
4069}
[860]4070
[1580]4071}
Note: See TracBrowser for help on using the repository browser.