source: GTP/trunk/Lib/Vis/Preprocessing/src/VspKdTree.cpp @ 744

Revision 744, 59.5 KB checked in by mattausch, 18 years ago (diff)
RevLine 
[408]1
2// ================================================================
3// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4// ****************************************************************
5/**
6   The KD tree based LSDS
7*/
8// Initial coding by
9/**
10   @author Jiri Bittner
11*/
12
13// Standard headers
[406]14#include <stack>
[408]15#include <queue>
[406]16#include <algorithm>
[408]17#include <fstream>
18#include <string>
19
20#include "VspKdTree.h"
21
[406]22#include "Environment.h"
[408]23#include "VssRay.h"
24#include "Intersectable.h"
25#include "Ray.h"
[420]26#include "RayInfo.h"
[453]27#include "ViewCellsManager.h"
[462]28#include "ViewCellBsp.h"
[406]29
[418]30// Static variables
[462]31int VspKdLeaf::sMailId = 0;
[580]32int VspKdMergeCandidate::sMaxPvsSize = 150;
33float VspKdMergeCandidate::sOverallCost = Limits::Small;
[428]34#define USE_FIXEDPOINT_T 0
35
[420]36/// Adds object to the pvs of the front and back node
[409]37inline void AddObject2Pvs(Intersectable *object,
38                                                  const int side,
39                                                  int &pvsBack,
40                                                  int &pvsFront)
[406]41{
[408]42        if (!object)
43                return;
[465]44
45        if (side <= 0)
[409]46        {
[465]47                if (!object->Mailed() && !object->Mailed(2))
[409]48                {
49                        ++ pvsBack;
50
[408]51                        if (object->Mailed(1))
52                                object->Mail(2);
53                        else
54                                object->Mail();
55                }
56        }
[465]57
[428]58        if (side >= 0)
[409]59        {
[465]60                if (!object->Mailed(1) && !object->Mailed(2))
[409]61                {
62                        ++ pvsFront;
[428]63
[408]64                        if (object->Mailed())
65                                object->Mail(2);
66                        else
67                                object->Mail(1);
68                }
69        }
[406]70}
71
[418]72/**************************************************************/
[462]73/*                class VspKdNode implementation          */
[420]74/**************************************************************/
75
76// Inline constructor
[500]77VspKdNode::VspKdNode(VspKdNode *p):
[465]78mParent(p), mAxis(-1), mDepth(p ? p->mDepth + 1 : 0)
[420]79{}
80
[465]81VspKdNode::~VspKdNode()
[480]82{
83};
[420]84
[500]85inline VspKdNode *VspKdNode::GetParent() const
[420]86{
87        return mParent;
88}
89
[500]90inline void VspKdNode::SetParent(VspKdNode *p)
[420]91{
92        mParent = p;
93}
94
[465]95bool VspKdNode::IsLeaf() const
96{
97        return mAxis == -1;
[420]98}
[465]99
100int VspKdNode::GetAccessTime()
[420]101{
102        return 0x7FFFFFF;
103}
104
105/**************************************************************/
[462]106/*                VspKdInterior implementation            */
[418]107/**************************************************************/
[408]108
[462]109VspKdInterior::VspKdInterior(VspKdInterior *p):
110VspKdNode(p), mBack(NULL), mFront(NULL), mAccesses(0), mLastAccessTime(-1)
[418]111{
112}
113
[462]114int VspKdInterior::GetAccessTime()
[418]115{
[419]116        return mLastAccessTime;
[418]117}
118
[465]119void VspKdInterior::SetupChildLinks(VspKdNode *b, VspKdNode *f)
[418]120{
121        mBack = b;
122        mFront = f;
[420]123        b->SetParent(this);
124        f->SetParent(this);
[418]125}
126
[465]127void VspKdInterior::ReplaceChildLink(VspKdNode *oldChild,
[500]128                                                                         VspKdNode *newChild)
[418]129{
130        if (mBack == oldChild)
131                mBack = newChild;
132        else
133                mFront = newChild;
134}
135
[465]136int VspKdInterior::Type() const
137{
138        return EInterior;
[418]139}
[465]140
[462]141VspKdInterior::~VspKdInterior()
[418]142{
143        DEL_PTR(mBack);
144        DEL_PTR(mFront);
145}
[465]146
147void VspKdInterior::Print(ostream &s) const
[418]148{
[423]149        switch (mAxis)
150        {
151        case 0:
152                s << "x "; break;
153        case 1:
154                s << "y "; break;
155        case 2:
156                s << "z "; break;
157        }
158
[419]159        s << mPosition << " ";
[418]160
161        mBack->Print(s);
162        mFront->Print(s);
163}
[465]164
165int VspKdInterior::ComputeRayIntersection(const RayInfo &rayData, float &t)
[418]166{
[419]167        return rayData.ComputeRayIntersection(mAxis, mPosition, t);
[418]168}
169
[462]170VspKdNode *VspKdInterior::GetBack() const
[419]171{
172        return mBack;
173}
[418]174
[462]175VspKdNode *VspKdInterior::GetFront() const
[406]176{
[419]177        return mFront;
178}
[406]179
[420]180
[486]181/*******************************************************************/
182/*                   class VspKdLeaf implementation                */
183/*******************************************************************/
[462]184
185
[500]186VspKdLeaf::VspKdLeaf(VspKdNode *p,
[472]187                                         const int nRays,
188                                         const int maxCostMisses):
189VspKdNode(p), mRays(), mPvsSize(0), mValidPvs(false), mViewCell(NULL),
[495]190mMaxCostMisses(maxCostMisses), mPvs(NULL), mVolume(0)
[420]191{
192        mRays.reserve(nRays);
193}
194
[465]195VspKdLeaf::~VspKdLeaf()
[462]196{
[495]197        DEL_PTR(mPvs);
[462]198}
[420]199
[472]200
201int VspKdLeaf::GetMaxCostMisses()
202{
203        return mMaxCostMisses;
204}
205
206
[465]207int VspKdLeaf::Type() const
208{
209        return ELeaf;
[420]210}
211
[465]212void VspKdLeaf::Print(ostream &s) const
[420]213{
214        s << endl << "L: r = " << (int)mRays.size() << endl;
215};
216
[465]217void VspKdLeaf::AddRay(const RayInfo &data)
[420]218{
219        mValidPvs = false;
220        mRays.push_back(data);
221        data.mRay->Ref();
222}
223
[465]224int VspKdLeaf::GetPvsSize() const
[420]225{
226        return mPvsSize;
227}
228
[465]229void VspKdLeaf::SetPvsSize(const int s)
[420]230{
231        mPvsSize = s;
232}
233
[462]234void VspKdLeaf::Mail()
[465]235{
236        mMailbox = sMailId;
[420]237}
238
[465]239void VspKdLeaf::NewMail()
240{
241        ++ sMailId;
[420]242}
243
[465]244bool VspKdLeaf::Mailed() const
245{
246        return mMailbox == sMailId;
[420]247}
248
[462]249bool VspKdLeaf::Mailed(const int mail) const
[420]250{
[462]251        return mMailbox == mail;
[420]252}
253
[462]254int VspKdLeaf::GetMailbox() const
[420]255{
[462]256        return mMailbox;
257}
258
[465]259float VspKdLeaf::GetAvgRayContribution() const
[462]260{
[420]261        return GetPvsSize() / ((float)mRays.size() + Limits::Small);
262}
263
[465]264float VspKdLeaf::GetSqrRayContribution() const
[420]265{
266        return sqr(GetAvgRayContribution());
267}
268
[462]269RayInfoContainer &VspKdLeaf::GetRays()
[426]270{
271        return mRays;
272}
[428]273
[462]274void VspKdLeaf::ExtractPvs(ObjectContainer &objects) const
[428]275{
276        RayInfoContainer::const_iterator it, it_end = mRays.end();
277
278        for (it = mRays.begin(); it != it_end; ++ it)
279        {
280                if ((*it).mRay->mTerminationObject)
281                        objects.push_back((*it).mRay->mTerminationObject);
282                if ((*it).mRay->mOriginObject)
283                        objects.push_back((*it).mRay->mOriginObject);
284        }
285}
286
[462]287void VspKdLeaf::GetRays(VssRayContainer &rays)
[428]288{
289        RayInfoContainer::const_iterator it, it_end = mRays.end();
290
291        for (it = mRays.begin(); it != mRays.end(); ++ it)
292                rays.push_back((*it).mRay);
293}
294
[462]295VspKdViewCell *VspKdLeaf::GetViewCell()
[453]296{
297        return mViewCell;
298}
299
[462]300void VspKdLeaf::SetViewCell(VspKdViewCell *viewCell)
301{
302        mViewCell = viewCell;
303}
304
305
306void VspKdLeaf::UpdatePvsSize()
307{
[465]308        if (!mValidPvs)
[462]309        {
310                Intersectable::NewMail();
311                int pvsSize = 0;
[465]312                for(RayInfoContainer::iterator ri = mRays.begin();
[462]313                        ri != mRays.end(); ++ ri)
314                {
[465]315                        if ((*ri).mRay->IsActive())
[462]316                        {
317                                Intersectable *object;
318#if BIDIRECTIONAL_RAY
319                                object = (*ri).mRay->mOriginObject;
[465]320
321                                if (object && !object->Mailed())
[462]322                                {
323                                        ++ pvsSize;
324                                        object->Mail();
325                                }
326#endif
327                                object = (*ri).mRay->mTerminationObject;
[465]328                                if (object && !object->Mailed())
[462]329                                {
330                                        ++ pvsSize;
331                                        object->Mail();
332                                }
333                        }
334                }
335
336                mPvsSize = pvsSize;
337                mValidPvs = true;
338        }
339}
340
[420]341
[500]342/***************************************************************/
343/*          class VspKdIntermediate implementation             */
344/***************************************************************/
[420]345
[500]346
347VspKdIntermediate::VspKdIntermediate(VspKdInterior *p):
348VspKdNode(p), mBack(NULL), mFront(NULL), mAccesses(0), mLastAccessTime(-1)
349{
350}
351
352
353VspKdIntermediate::~VspKdIntermediate()
354{
355        DEL_PTR(mFront);
356        DEL_PTR(mBack);
357}
358   
359       
360int VspKdIntermediate::Type() const
361{
362        return EIntermediate;
363}
364
365
366VspKdLeaf *VspKdIntermediate::GetBack() const
367{
368        return mBack;
369}
370
371
372VspKdLeaf *VspKdIntermediate::GetFront() const
373{
374        return mFront;
375}
376
377
378void VspKdIntermediate::Print(ostream &s) const
379{
380        s << endl << mSplitPlane << endl;
381}
382
383
384void VspKdIntermediate::SetupChildLinks(VspKdLeaf *b, VspKdLeaf *f)
385{
386        mBack = b;
387        mFront = f;
388
389        b->SetParent(this);
390        f->SetParent(this);
391}
392
393
394int VspKdIntermediate::ComputeRayIntersection(const RayInfo &rayData, float &t)
395{
396        return rayData.ComputeRayIntersection(mSplitPlane, t);
397}
398
399
400/*************************************************************/
401/*              class VspKdTree implementation               */
402/*************************************************************/
403
404
[482]405VspKdTree::VspKdTree(): mOnlyDrivingAxis(false)
[465]406{
[421]407        environment->GetIntValue("VspKdTree.Termination.maxDepth", mTermMaxDepth);
408        environment->GetIntValue("VspKdTree.Termination.minPvs", mTermMinPvs);
409        environment->GetIntValue("VspKdTree.Termination.minRays", mTermMinRays);
410        environment->GetFloatValue("VspKdTree.Termination.maxRayContribution", mTermMaxRayContribution);
411        environment->GetFloatValue("VspKdTree.Termination.maxCostRatio", mTermMaxCostRatio);
[422]412        environment->GetFloatValue("VspKdTree.Termination.minSize", mTermMinSize);
[421]413
[472]414        environment->GetIntValue("VspKdTree.Termination.missTolerance", mTermMissTolerance);
415
[421]416        mTermMinSize = sqr(mTermMinSize);
417
[422]418        environment->GetFloatValue("VspKdTree.epsilon", mEpsilon);
419        environment->GetFloatValue("VspKdTree.ct_div_ci", mCtDivCi);
[465]420
[422]421        environment->GetFloatValue("VspKdTree.maxTotalMemory", mMaxTotalMemory);
422        environment->GetFloatValue("VspKdTree.maxStaticMemory", mMaxStaticMemory);
[465]423
[422]424        environment->GetIntValue("VspKdTree.accessTimeThreshold", mAccessTimeThreshold);
425        environment->GetIntValue("VspKdTree.minCollapseDepth", mMinCollapseDepth);
[408]426
[465]427        //environment->GetIntValue("ViewCells.maxViewCells", mMaxViewCells);
[462]428
[472]429        environment->GetIntValue("VspKdTree.PostProcess.minViewCells", mMergeMinViewCells);
430        environment->GetFloatValue("VspKdTree.PostProcess.maxCostRatio", mMergeMaxCostRatio);
[465]431        environment->GetIntValue("VspKdTree.PostProcess.maxPvsSize",
[580]432                                                         VspKdMergeCandidate::sMaxPvsSize);
[462]433
[485]434        environment->GetBoolValue("VspKdTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
[587]435        environment->GetIntValue("VspKdTree.Termination.maxViewCells", mMaxViewCells);
[465]436
[485]437        //-- output
[423]438        Debug << "======= vsp kd tree options ========" << endl;
439        Debug << "max depth: "<< mTermMaxDepth << endl;
440        Debug << "min pvs: "<< mTermMinPvs << endl;
441        Debug << "min rays: "<< mTermMinRays << endl;
442        Debug << "max ray contribution: "<< mTermMaxRayContribution << endl;
[473]443        Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
444        Debug << "min size: " << mTermMinSize << endl;
[423]445
[485]446        Debug << "split type: ";
447
448        //-- split type
449        char sname[128];
450        environment->GetStringValue("VspKdTree.splitType", sname);
451        string name(sname);
452
[409]453        if (name.compare("regular") == 0)
[423]454        {
[485]455                Debug << "regular split" << endl;
[409]456                splitType = ESplitRegular;
[423]457        }
[409]458        else
459        {
460                if (name.compare("heuristic") == 0)
[423]461                {
[485]462                        Debug << "heuristic split" << endl;
[409]463                        splitType = ESplitHeuristic;
[423]464                }
[465]465                else
[409]466                {
467                        cerr << "Invalid VspKdTree split type " << name << endl;
[408]468                        exit(1);
469                }
[409]470        }
[408]471
[418]472        mRoot = NULL;
[422]473        mSplitCandidates = new vector<SortableEntry>;
[406]474}
475
[408]476
477VspKdTree::~VspKdTree()
[406]478{
[419]479        DEL_PTR(mRoot);
[449]480        DEL_PTR(mSplitCandidates);
[408]481}
[406]482
[418]483void VspKdStatistics::Print(ostream &app) const
[408]484{
[426]485        app << "===== VspKdTree statistics ===============\n";
[408]486
[426]487        app << "#N_RAYS ( Number of rays )\n"
488                << rays << endl;
[465]489
[426]490        app << "#N_INITPVS ( Initial PVS size )\n"
491                << initialPvsSize << endl;
[465]492
[426]493        app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
[406]494
[426]495        app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
[406]496
[667]497        app << "#N_SPLITS ( Number of splits in axes x y z)\n";
[465]498
[482]499        for (int i = 0; i < 3; ++ i)
500                app << splits[i] << " ";
[426]501        app << endl;
[406]502
[426]503        app << "#N_RAYREFS ( Number of rayRefs )\n" <<
504                rayRefs << "\n";
[408]505
[426]506        app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
507                rayRefs / (double)rays << "\n";
[408]508
[426]509        app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
510                rayRefs / (double)Leaves() << "\n";
[408]511
[426]512        app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
513                maxRayRefs << "\n";
[408]514
515  //  app << setprecision(4);
516
[426]517        app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n" <<
518                maxDepthNodes * 100 / (double)Leaves() << endl;
[408]519
[426]520        app << "#N_PMINPVSLEAVES  ( Percentage of leaves with mininimal PVS )\n" <<
521                minPvsNodes * 100 / (double)Leaves() << endl;
[408]522
[426]523        app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minimal number of rays)\n" <<
524                minRaysNodes * 100 / (double)Leaves() << endl;
525
[408]526        app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
[426]527                minSizeNodes * 100 / (double)Leaves() << endl;
[408]528
[426]529        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n" <<
530                maxRayContribNodes * 100 / (double)Leaves() << endl;
[408]531
[482]532        app << "#N_PMAXRCOSTLEAVES  ( Percentage of leaves with maximal cost ratio )\n" <<
533                maxCostNodes * 100 / (double)Leaves() << endl;
534
[426]535        app << "#N_ADDED_RAYREFS  ( Number of dynamically added ray references )\n"<<
536                addedRayRefs << endl;
[408]537
[426]538        app << "#N_REMOVED_RAYREFS  ( Number of dynamically removed ray references )\n"<<
539                removedRayRefs << endl;
[408]540
[426]541        //  app << setprecision(4);
[408]542
[462]543        app << "#N_( Maximal PVS size / leaf)\n"
[426]544                << maxPvsSize << endl;
[408]545
[482]546        app << "#N_( Average PVS size / leaf)\n"
547                << (double) accPvsSize / (double)Leaves() << endl;
548
[426]549        app << "#N_CTIME  ( Construction time [s] )\n"
550                << Time() << " \n";
[408]551
[426]552        app << "===== END OF VspKdTree statistics ==========\n";
[406]553}
554
[408]555
[463]556void VspKdTree::CollectViewCells(ViewCellContainer &viewCells) const
557{
558        stack<VspKdNode *> nodeStack;
559        nodeStack.push(mRoot);
560
561        ViewCell::NewMail();
562
[465]563        while (!nodeStack.empty())
[463]564        {
565                VspKdNode *node = nodeStack.top();
566                nodeStack.pop();
567
[465]568                if (node->IsLeaf())
[463]569                {
570                        ViewCell *viewCell = dynamic_cast<VspKdLeaf *>(node)->mViewCell;
571
[465]572                        if (!viewCell->Mailed())
[463]573                        {
574                                viewCell->Mail();
575                                viewCells.push_back(viewCell);
576                        }
577                }
[465]578                else
[463]579                {
580                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
581
582                        nodeStack.push(interior->GetFront());
583                        nodeStack.push(interior->GetBack());
584                }
585        }
[408]586}
[406]587
[462]588void VspKdTree::Construct(const VssRayContainer &rays,
589                                                  AxisAlignedBox3 *forcedBoundingBox)
[452]590{
[414]591        mStat.Start();
[465]592
[422]593        mMaxMemory = mMaxStaticMemory;
[419]594        DEL_PTR(mRoot);
[408]595
[462]596        mRoot = new VspKdLeaf(NULL, (int)rays.size());
[408]597
[423]598        // first construct a leaf that will get subdivided
[462]599        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(mRoot);
[465]600
[414]601        mStat.nodes = 1;
[448]602        mBox.Initialize();
[452]603
[448]604        //-- compute bounding box
[409]605        if (forcedBoundingBox)
[419]606                mBox = *forcedBoundingBox;
[448]607        else
608                for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
609                {
[452]610                        if ((*ri)->mOriginObject)
611                mBox.Include((*ri)->GetOrigin());
612                        if ((*ri)->mTerminationObject)
[465]613                                mBox.Include((*ri)->GetTermination());
[448]614                }
[408]615
[422]616        Debug << "bbox = " << mBox << endl;
[448]617
618        for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
619        {
620                //leaf->AddRay(RayInfo(*ri));
[463]621                VssRay *ray = *ri;
622                float minT, maxT;
[482]623               
624                // TODO: not very efficient to implictly cast between rays types!
[463]625                if (mBox.GetRaySegment(*ray, minT, maxT))
626                {
627                        float len = ray->Length();
[465]628
629                        if (!len)
[463]630                                len = Limits::Small;
[465]631
632                        leaf->AddRay(RayInfo(ray, minT / len, maxT / len));
[448]633                }
634        }
[465]635
[420]636        mStat.rays = (int)leaf->mRays.size();
[408]637        leaf->UpdatePvsSize();
[465]638
[414]639        mStat.initialPvsSize = leaf->GetPvsSize();
[465]640
[409]641        // Subdivide();
[419]642        mRoot = Subdivide(TraversalData(leaf, mBox, 0));
[408]643
[465]644        if (mSplitCandidates)
[409]645        {
[422]646                // force release of this vector
647                delete mSplitCandidates;
648                mSplitCandidates = new vector<SortableEntry>;
[409]649        }
[465]650
[414]651        mStat.Stop();
[408]652
[414]653        mStat.Print(cout);
[423]654        Debug << "#Total memory=" << GetMemUsage() << endl;
[408]655}
656
657
658
[462]659VspKdNode *VspKdTree::Subdivide(const TraversalData &tdata)
[408]660{
[462]661        VspKdNode *result = NULL;
[408]662
[409]663        priority_queue<TraversalData> tStack;
[423]664        //stack<TraversalData> tStack;
[465]665
[409]666        tStack.push(tdata);
[406]667
[409]668        AxisAlignedBox3 backBox;
669        AxisAlignedBox3 frontBox;
[465]670
[408]671        int lastMem = 0;
[423]672
[465]673        while (!tStack.empty())
[409]674        {
[408]675                float mem = GetMemUsage();
[465]676
677                if (lastMem / 10 != ((int)mem) / 10)
[409]678                {
[425]679                        Debug << mem << " MB" << endl;
[408]680                }
681                lastMem = (int)mem;
[465]682
683                if (1 && mem > mMaxMemory)
[409]684                {
[423]685                        Debug << "memory limit reached: " << mem << endl;
[409]686                        // count statistics on unprocessed leafs
[465]687                        while (!tStack.empty())
[409]688                        {
[408]689                                EvaluateLeafStats(tStack.top());
690                                tStack.pop();
[409]691                        }
692                        break;
693                }
[465]694
[409]695                TraversalData data = tStack.top();
[465]696                tStack.pop();
697
[472]698                VspKdNode *node = SubdivideNode(dynamic_cast<VspKdLeaf *>(data.mNode),
699                                                                                data.mBox, backBox,     frontBox);
[409]700                if (result == NULL)
701                        result = node;
[465]702
703                if (!node->IsLeaf())
[409]704                {
[462]705                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
[408]706
[409]707                        // push the children on the stack
[419]708                        tStack.push(TraversalData(interior->GetBack(), backBox, data.mDepth + 1));
709                        tStack.push(TraversalData(interior->GetFront(), frontBox, data.mDepth + 1));
[465]710                }
711                else
[409]712                {
713                        EvaluateLeafStats(data);
714                }
715        }
716
717        return result;
[408]718}
[406]719
[408]720
721// returns selected plane for subdivision
[462]722int VspKdTree::SelectPlane(VspKdLeaf *leaf,
[424]723                                                   const AxisAlignedBox3 &box,
724                                                   float &position,
725                                                   int &raysBack,
726                                                   int &raysFront,
727                                                   int &pvsBack,
728                                                   int &pvsFront)
[408]729{
[428]730        int axis = 0;
731        float costRatio = 0;
[465]732
733        if (splitType == ESplitRegular)
[422]734        {
[408]735                costRatio = BestCostRatioRegular(leaf,
[480]736                                                                                 box,
[410]737                                                                                 axis,
738                                                                                 position,
739                                                                                 raysBack,
740                                                                                 raysFront,
741                                                                                 pvsBack,
742                                                                                 pvsFront);
[465]743        }
[423]744        else if (splitType == ESplitHeuristic)
[410]745        {
[472]746                costRatio = BestCostRatioHeuristic(leaf,
[480]747                                                                                   box,
[472]748                                                                                   axis,
749                                                                                   position,
750                                                                                   raysBack,
751                                                                                   raysFront,
752                                                                                   pvsBack,
753                                                                                   pvsFront);
[423]754        }
[465]755        else
[423]756        {
757                cerr << "VspKdTree: Unknown split heuristics\n";
758                exit(1);
759        }
[408]760
[465]761        if (costRatio > mTermMaxCostRatio)
[410]762        {
[482]763                //Debug << "Too big cost ratio " << costRatio << endl;
764                ++ leaf->mMaxCostMisses;
765                if (leaf->mMaxCostMisses > mTermMissTolerance)
766                {
767                        ++ mStat.maxCostNodes;
768                        return -1;
769                }
[408]770        }
771
[428]772        if (0)
[465]773                Debug << "pvs=" << leaf->mPvsSize
774                          << " rays=" << (int)leaf->mRays.size()
775                          << " rc=" << leaf->GetAvgRayContribution()
[428]776                          << " axis=" << axis << endl;
[465]777
[408]778        return axis;
[406]779}
780
781
[465]782
[462]783float VspKdTree::EvalCostRatio(VspKdLeaf *leaf,
[480]784                                                           const AxisAlignedBox3 &box,
[425]785                                                           const int axis,
786                                                           const float position,
787                                                           int &raysBack,
788                                                           int &raysFront,
789                                                           int &pvsBack,
790                                                           int &pvsFront)
[406]791{
[408]792        raysBack = 0;
793        raysFront = 0;
794        pvsFront = 0;
795        pvsBack = 0;
796
797        Intersectable::NewMail(3);
[465]798
[408]799        // eval pvs size
[482]800        const int pvsSize = leaf->GetPvsSize();
[408]801
[419]802        // this is the main ray classification loop!
[480]803        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[463]804                        ri != leaf->mRays.end(); ++ ri)
[480]805        {
806                if (!(*ri).mRay->IsActive())
807                        continue;
[465]808
[480]809                // determine the side of this ray with respect to the plane
810                float t;
811                int side = (*ri).ComputeRayIntersection(axis, position, t);
812                       
[419]813                if (side <= 0)
814                        ++ raysBack;
[465]815
[419]816                if (side >= 0)
817                        ++ raysFront;
[465]818
[419]819                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
[465]820        }
[408]821
[482]822        //-- only one of these options should be one
823
824        if (0) //-- only pvs
[425]825        {
[428]826                const float sum = float(pvsBack + pvsFront);
[482]827                const float oldCost = (float)pvsSize;
[425]828
829                return sum / oldCost;
830        }
[482]831
832        //-- pvs + probability heuristics
833        float pBack, pFront, pOverall;
834
[542]835        if (1)
[428]836        {
[482]837                // box length substitute for probability
838                const float minBox = box.Min(axis);
839                const float maxBox = box.Max(axis);
[480]840
[482]841                pBack = position - minBox;
842                pFront = maxBox - position;
843                pOverall = maxBox - minBox;
844        }
[465]845
[542]846        if (0) //-- area substitute for probability
[482]847        {
[495]848                pOverall = box.SurfaceArea();
849
[482]850                const bool useMidSplit = true;
[495]851                //const bool useMidSplit = false;       
852                                       
[482]853                if (!useMidSplit)
854                {
855                        Vector3 pos = box.Max(); pos[axis] = position;
856                        pBack = AxisAlignedBox3(box.Min(), pos).SurfaceArea();
[480]857
[482]858                        pos = box.Min(); pos[axis] = position;
859                        pFront = AxisAlignedBox3(pos, box.Max()).SurfaceArea();
860                }
861                else
862                {
863                        //-- simplified computation for mid split
864                        const int axis2 = (axis + 1) % 3;
865                        const int axis3 = (axis + 2) % 3;
[480]866
[482]867                        const float faceArea =
868                                (box.Max(axis2) - box.Min(axis2)) *
869                                (box.Max(axis3) - box.Min(axis3));
[408]870
[482]871                        pBack = pFront = pOverall * 0.5f + faceArea;
872                }
873        }
[465]874
[482]875        //-- ray density substitute for probability
[483]876        if (0)
[482]877        {
878                pBack = (float)raysBack;
879                pFront = (float)raysFront;
880                pOverall = (float)leaf->mRays.size();
881        }
[465]882
[482]883        //Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
884        //Debug << pFront << " " << pBack << " " << pOverall << endl;
885
886        // float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
887        const float newCost = pvsBack * pBack + pvsFront * pFront;
888        //  float oldCost = leaf->mRays.size();
889        const float oldCost = (float)pvsSize * pOverall;
890
891        return  (mCtDivCi + newCost) / oldCost;
[406]892}
893
[491]894
[462]895float VspKdTree::BestCostRatioRegular(VspKdLeaf *leaf,
[480]896                                                                          const AxisAlignedBox3 &box,
[422]897                                                                          int &axis,
898                                                                          float &position,
899                                                                          int &raysBack,
900                                                                          int &raysFront,
901                                                                          int &pvsBack,
902                                                                          int &pvsFront)
[408]903{
[422]904        int nRaysBack[3], nRaysFront[3];
905        int nPvsBack[3], nPvsFront[3];
906
907        float nPosition[3];
908        float nCostRatio[3];
[408]909        int bestAxis = -1;
[465]910
[480]911        //AxisAlignedBox3 sBox = GetBBox(leaf);
912        const int sAxis = box.Size().DrivingAxis();
[408]913
[465]914        for (axis = 0; axis < 3; ++ axis)
[410]915        {
[465]916                if (!mOnlyDrivingAxis || axis == sAxis)
[410]917                {
[465]918
[480]919                        nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
[465]920
[408]921                        nCostRatio[axis] = EvalCostRatio(leaf,
[480]922                                                                                         box,
[410]923                                                                                         axis,
924                                                                                         nPosition[axis],
925                                                                                         nRaysBack[axis],
926                                                                                         nRaysFront[axis],
927                                                                                         nPvsBack[axis],
928                                                                                         nPvsFront[axis]);
929                        if (bestAxis == -1)
[452]930                        {
[408]931                                bestAxis = axis;
[452]932                        }
933                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
934                        {
[465]935                                /*Debug << "pvs front " << nPvsBack[axis]
[452]936                                          << " pvs back " << nPvsFront[axis]
937                                          << " overall pvs " << leaf->GetPvsSize() << endl;*/
938                                bestAxis = axis;
939                        }
[465]940
[408]941                }
942        }
943
[428]944        //-- assign best axis
[408]945        axis = bestAxis;
946        position = nPosition[bestAxis];
947
948        raysBack = nRaysBack[bestAxis];
949        raysFront = nRaysFront[bestAxis];
950
951        pvsBack = nPvsBack[bestAxis];
952        pvsFront = nPvsFront[bestAxis];
[465]953
[408]954        return nCostRatio[bestAxis];
955}
956
[462]957float VspKdTree::BestCostRatioHeuristic(VspKdLeaf *leaf,
[480]958                                                                                const AxisAlignedBox3 &box,
[410]959                                                                                int &axis,
960                                                                                float &position,
961                                                                                int &raysBack,
962                                                                                int &raysFront,
963                                                                                int &pvsBack,
964                                                                                int &pvsFront)
[406]965{
[480]966        //AxisAlignedBox3 box = GetBBox(leaf);
[408]967        //      AxisAlignedBox3 dirBox = GetDirBBox(node);
[465]968
[408]969        axis = box.Size().DrivingAxis();
[465]970
[408]971        SortSplitCandidates(leaf, axis);
[465]972
[408]973        // go through the lists, count the number of objects left and right
974        // and evaluate the following cost funcion:
975        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
[465]976
[420]977        int rl = 0, rr = (int)leaf->mRays.size();
978        int pl = 0, pr = leaf->GetPvsSize();
[424]979
[408]980        float minBox = box.Min(axis);
981        float maxBox = box.Max(axis);
982        float sizeBox = maxBox - minBox;
[465]983
[419]984        float minBand = minBox + 0.1f*(maxBox - minBox);
985        float maxBand = minBox + 0.9f*(maxBox - minBox);
[465]986
[408]987        float sum = rr*sizeBox;
[411]988        float minSum = 1e20f;
[465]989
[408]990        Intersectable::NewMail();
[424]991
[408]992        // set all object as belonging to the fron pvs
[420]993        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[424]994                ri != leaf->mRays.end(); ++ ri)
[410]995        {
996                if ((*ri).mRay->IsActive())
997                {
[408]998                        Intersectable *object = (*ri).mRay->mTerminationObject;
[465]999
[408]1000                        if (object)
[465]1001                                if (!object->Mailed())
[410]1002                                {
[408]1003                                        object->Mail();
1004                                        object->mCounter = 1;
[465]1005                                }
[410]1006                                else
1007                                        ++ object->mCounter;
[408]1008                }
[410]1009        }
[465]1010
[408]1011        Intersectable::NewMail();
[465]1012
[422]1013        for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
[465]1014                ci < mSplitCandidates->end(); ++ ci)
[410]1015        {
[408]1016                VssRay *ray;
[465]1017
1018                switch ((*ci).type)
[410]1019                {
[465]1020                        case SortableEntry::ERayMin:
[410]1021                                {
1022                                        ++ rl;
[482]1023                                        ray = (*ci).ray;
[410]1024                                        Intersectable *object = ray->mTerminationObject;
[465]1025                                        if (object && !object->Mailed())
[410]1026                                        {
1027                                                object->Mail();
1028                                                ++ pl;
1029                                        }
1030                                        break;
1031                                }
[465]1032                        case SortableEntry::ERayMax:
[410]1033                                {
1034                                        -- rr;
[465]1035
[482]1036                                        ray = (*ci).ray;
[410]1037                                        Intersectable *object = ray->mTerminationObject;
[465]1038
[410]1039                                        if (object)
1040                                        {
1041                                                if (-- object->mCounter == 0)
1042                                                -- pr;
1043                                        }
[465]1044
[410]1045                                        break;
1046                                }
[408]1047                }
[465]1048
1049                if ((*ci).value > minBand && (*ci).value < maxBand)
[410]1050                {
[408]1051                        sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);
[465]1052
[410]1053                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
1054                        // cout<<"cost= "<<sum<<endl;
[465]1055
1056                        if (sum < minSum)
[410]1057                        {
[408]1058                                minSum = sum;
1059                                position = (*ci).value;
[465]1060
[408]1061                                raysBack = rl;
1062                                raysFront = rr;
[465]1063
[408]1064                                pvsBack = pl;
1065                                pvsFront = pr;
[465]1066
[410]1067                        }
1068                }
1069        }
1070
[419]1071        float oldCost = (float)leaf->GetPvsSize();
[422]1072        float newCost = mCtDivCi + minSum / sizeBox;
[410]1073        float ratio = newCost / oldCost;
[465]1074
[422]1075        //Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1076        //     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
[465]1077
[408]1078        return ratio;
1079}
1080
[462]1081void VspKdTree::SortSplitCandidates(VspKdLeaf *node,
[410]1082                                                                        const int axis)
[408]1083{
[422]1084        mSplitCandidates->clear();
[465]1085
[420]1086        int requestedSize = 2 * (int)(node->mRays.size());
[410]1087        // creates a sorted split candidates array
[422]1088        if (mSplitCandidates->capacity() > 500000 &&
[465]1089                requestedSize < (int)(mSplitCandidates->capacity()/10) )
[410]1090        {
[449]1091        DEL_PTR(mSplitCandidates);
[422]1092                mSplitCandidates = new vector<SortableEntry>;
[410]1093        }
[465]1094
[422]1095        mSplitCandidates->reserve(requestedSize);
[408]1096
[465]1097        // insert all queries
[420]1098        for(RayInfoContainer::const_iterator ri = node->mRays.begin();
[465]1099                ri < node->mRays.end(); ++ ri)
[410]1100        {
1101                bool positive = (*ri).mRay->HasPosDir(axis);
[465]1102                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
[482]1103                                                                                                  (*ri).ExtrapOrigin(axis), (*ri).mRay));
[465]1104
[422]1105                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
[482]1106                                                                                                  (*ri).ExtrapTermination(axis), (*ri).mRay));
[410]1107        }
[465]1108
[422]1109        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
[406]1110}
1111
[408]1112
[410]1113void VspKdTree::EvaluateLeafStats(const TraversalData &data)
[406]1114{
[410]1115        // the node became a leaf -> evaluate stats for leafs
[462]1116        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(data.mNode);
[408]1117
[426]1118        if (leaf->GetPvsSize() > mStat.maxPvsSize)
1119                mStat.maxPvsSize = leaf->GetPvsSize();
1120
1121        if ((int)(leaf->mRays.size()) > mStat.maxRayRefs)
1122                mStat.maxRayRefs = (int)leaf->mRays.size();
1123
[482]1124        mStat.rays += (int)leaf->mRays.size();
[426]1125
[421]1126        if (data.mDepth >= mTermMaxDepth)
[744]1127        {
[413]1128                ++ mStat.maxDepthNodes;
[744]1129        }
[465]1130
[421]1131        if (leaf->GetPvsSize() < mTermMinPvs)
[413]1132                ++ mStat.minPvsNodes;
[408]1133
[482]1134        mStat.accPvsSize += leaf->GetPvsSize();
1135
[426]1136        if ((int)leaf->GetRays().size() < mTermMinRays)
[413]1137                ++ mStat.minRaysNodes;
[408]1138
[426]1139        if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
[413]1140                ++ mStat.maxRayContribNodes;
[465]1141
1142        if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
[413]1143                ++ mStat.minSizeNodes;
[408]1144}
1145
1146
[465]1147inline bool VspKdTree::TerminationCriteriaMet(const VspKdLeaf *leaf,
[421]1148                                                                                          const AxisAlignedBox3 &box) const
1149{
[465]1150        return ((leaf->GetPvsSize() < mTermMinPvs) ||
[426]1151                    (leaf->mRays.size() < mTermMinRays) ||
[482]1152                        (leaf->GetAvgRayContribution() > mTermMaxRayContribution ) ||
[465]1153                        (leaf->mDepth >= mTermMaxDepth) ||
[587]1154                        (mStat.Leaves() >= mMaxViewCells) ||
[423]1155                        (SqrMagnitude(box.Size()) <= mTermMinSize));
[421]1156}
[408]1157
[421]1158
[462]1159VspKdNode *VspKdTree::SubdivideNode(VspKdLeaf *leaf,
[472]1160                                                                        const AxisAlignedBox3 &box,
1161                                                                        AxisAlignedBox3 &backBBox,
1162                                                                        AxisAlignedBox3 &frontBBox)
[408]1163{
[421]1164        if (TerminationCriteriaMet(leaf, box))
[410]1165        {
[465]1166                if (1 && leaf->mDepth >= mTermMaxDepth)
[422]1167                {
[428]1168                        Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
1169                        Debug << "Bbox: " << GetBBox(leaf) << endl;
[408]1170                }
[428]1171                //Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
[465]1172
[408]1173                return leaf;
1174        }
[465]1175
[408]1176        float position;
1177        // first count ray sides
[410]1178        int raysBack;
1179        int raysFront;
[408]1180        int pvsBack;
1181        int pvsFront;
[465]1182
[410]1183        // select subdivision axis
[425]1184        const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
[426]1185        //Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" <<       pvsFront << endl;
[408]1186
[465]1187        if (axis == -1)
[482]1188                return leaf;
1189       
[413]1190        mStat.nodes += 2;
[482]1191        ++ mStat.splits[axis];
[408]1192
[410]1193        // add the new nodes to the tree
[500]1194        VspKdInterior *node =
1195                new VspKdInterior(dynamic_cast<VspKdInterior *>(leaf->mParent));
[406]1196
[419]1197        node->mAxis = axis;
1198        node->mPosition = position;
1199        node->mBox = box;
[465]1200
[410]1201        backBBox = box;
1202        frontBBox = box;
[406]1203
[482]1204        VspKdLeaf *back = new VspKdLeaf(node, raysBack, leaf->GetMaxCostMisses());
[408]1205        back->SetPvsSize(pvsBack);
[482]1206        VspKdLeaf *front = new VspKdLeaf(node, raysFront, leaf->GetMaxCostMisses());
[408]1207        front->SetPvsSize(pvsFront);
[406]1208
[410]1209        // replace a link from node's parent
[419]1210        if (leaf->mParent)
[500]1211        {
1212                VspKdInterior *parent = dynamic_cast<VspKdInterior *>(leaf->mParent);
1213                parent->ReplaceChildLink(leaf, node);
1214        }
[501]1215
[410]1216        // and setup child links
1217        node->SetupChildLinks(back, front);
[465]1218
1219        if (axis <= VspKdNode::SPLIT_Z)
[410]1220        {
[408]1221                backBBox.SetMax(axis, position);
[410]1222                frontBBox.SetMin(axis, position);
[465]1223
[420]1224                for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1225                                ri != leaf->mRays.end(); ++ ri)
[410]1226                {
[465]1227                        if ((*ri).mRay->IsActive())
[410]1228                        {
[408]1229                                // first unref ray from the former leaf
1230                                (*ri).mRay->Unref();
[463]1231                                  float t;
[408]1232                                // determine the side of this ray with respect to the plane
[463]1233                                int side = node->ComputeRayIntersection(*ri, t);
[408]1234
[465]1235                                if (side == 0)
[410]1236                                {
[465]1237                                        if ((*ri).mRay->HasPosDir(axis))
[410]1238                                        {
[420]1239                                                back->AddRay(RayInfo((*ri).mRay,
1240                                                                                         (*ri).mMinT,
[463]1241                                                                                         t));
[420]1242                                                front->AddRay(RayInfo((*ri).mRay,
[463]1243                                                                                          t,
[420]1244                                                                                          (*ri).mMaxT));
[465]1245                                        }
1246                                        else
[410]1247                                        {
[420]1248                                                back->AddRay(RayInfo((*ri).mRay,
[463]1249                                                                                         t,
[420]1250                                                                                         (*ri).mMaxT));
1251                                                front->AddRay(RayInfo((*ri).mRay,
1252                                                                                          (*ri).mMinT,
[463]1253                                                                                          t));
[408]1254                                        }
[465]1255                                }
[410]1256                                else
[465]1257                                        if (side == 1)
[408]1258                                                front->AddRay(*ri);
1259                                        else
1260                                                back->AddRay(*ri);
[410]1261                        } else
[408]1262                                (*ri).mRay->Unref();
[410]1263                }
[465]1264        }
1265        else
[410]1266        {
1267                // rays front/back
[465]1268
[420]1269        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1270                        ri != leaf->mRays.end(); ++ ri)
[410]1271                {
[465]1272                        if ((*ri).mRay->IsActive())
[410]1273                        {
[408]1274                                // first unref ray from the former leaf
1275                                (*ri).mRay->Unref();
1276
1277                                int side;
1278                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1279                                        side = 1;
1280                                else
1281                                        side = -1;
[465]1282
[408]1283                                if (side == 1)
1284                                        front->AddRay(*ri);
1285                                else
[465]1286                                        back->AddRay(*ri);
1287                        }
[410]1288                        else
[408]1289                                (*ri).mRay->Unref();
[410]1290                }
1291        }
[465]1292
[410]1293        // update stats
[420]1294        mStat.rayRefs -= (int)leaf->mRays.size();
[414]1295        mStat.rayRefs += raysBack + raysFront;
[408]1296
[420]1297        DEL_PTR(leaf);
1298
[410]1299        return node;
[406]1300}
1301
[410]1302int VspKdTree::ReleaseMemory(const int time)
1303{
[462]1304        stack<VspKdNode *> tstack;
[406]1305
[410]1306        // find a node in the tree which subtree will be collapsed
[422]1307        int maxAccessTime = time - mAccessTimeThreshold;
[420]1308        int released = 0;
[406]1309
[419]1310        tstack.push(mRoot);
[406]1311
[423]1312        while (!tstack.empty())
[410]1313        {
[462]1314                VspKdNode *node = tstack.top();
[410]1315                tstack.pop();
[465]1316
1317                if (!node->IsLeaf())
[410]1318                {
[462]1319                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
[410]1320                        // cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
[465]1321                        if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
[410]1322                        {
1323                                released = CollapseSubtree(node, time);
1324                                break;
1325                        }
[465]1326                        if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
[410]1327                        {
[419]1328                                tstack.push(in->GetFront());
1329                                tstack.push(in->GetBack());
[465]1330                        }
1331                        else
[410]1332                        {
[419]1333                                tstack.push(in->GetBack());
1334                                tstack.push(in->GetFront());
[410]1335                        }
1336                }
1337        }
[406]1338
[465]1339        while (tstack.empty())
[410]1340        {
1341                // could find node to collaps...
1342                // cout<<"Could not find a node to release "<<endl;
1343                break;
1344        }
[465]1345
[410]1346        return released;
[408]1347}
[406]1348
1349
[462]1350VspKdNode *VspKdTree::SubdivideLeaf(VspKdLeaf *leaf,
[410]1351                                                                                const float sizeThreshold)
[408]1352{
[462]1353        VspKdNode *node = leaf;
[406]1354
[410]1355        AxisAlignedBox3 leafBBox = GetBBox(leaf);
[406]1356
[408]1357        static int pass = 0;
[410]1358        ++ pass;
[465]1359
[410]1360        // check if we should perform a dynamic subdivision of the leaf
[420]1361        if (// leaf->mRays.size() > (unsigned)termMinCost &&
[465]1362                (leaf->GetPvsSize() >= mTermMinPvs) &&
[410]1363                (SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
1364        {
[423]1365        // memory check and release
[465]1366                if (GetMemUsage() > mMaxTotalMemory)
[410]1367                        ReleaseMemory(pass);
[465]1368
[410]1369                AxisAlignedBox3 backBBox, frontBBox;
[406]1370
[410]1371                // subdivide the node
1372                node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
1373        }
1374        return node;
[406]1375}
1376
1377
1378
[423]1379void VspKdTree::UpdateRays(VssRayContainer &remove,
1380                                                   VssRayContainer &add)
[406]1381{
[462]1382        VspKdLeaf::NewMail();
[406]1383
[410]1384        // schedule rays for removal
1385        for(VssRayContainer::const_iterator ri = remove.begin();
[465]1386                ri != remove.end(); ++ ri)
[410]1387        {
[465]1388                (*ri)->ScheduleForRemoval();
[410]1389        }
[406]1390
[410]1391        int inactive = 0;
[406]1392
[465]1393        for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
[410]1394        {
1395                if ((*ri)->ScheduledForRemoval())
1396                        //      RemoveRay(*ri, NULL, false);
1397                        // !!! BUG - with true it does not work correctly - aggreated delete
1398                        RemoveRay(*ri, NULL, true);
1399                else
1400                        ++ inactive;
1401        }
[406]1402
[410]1403        //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
[465]1404        for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
[410]1405        {
1406                AddRay(*ri);
1407        }
[406]1408}
1409
1410
[410]1411void VspKdTree::RemoveRay(VssRay *ray,
[462]1412                                                  vector<VspKdLeaf *> *affectedLeaves,
[410]1413                                                  const bool removeAllScheduledRays)
[406]1414{
[410]1415        stack<RayTraversalData> tstack;
[406]1416
[420]1417        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[465]1418
[410]1419        RayTraversalData data;
[406]1420
[410]1421        // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
[465]1422        while (!tstack.empty())
[410]1423        {
1424                data = tstack.top();
1425                tstack.pop();
[406]1426
[465]1427                if (!data.mNode->IsLeaf())
[410]1428                {
1429                        // split the set of rays in two groups intersecting the
1430                        // two subtrees
1431                        TraverseInternalNode(data, tstack);
[465]1432        }
1433                else
[410]1434                {
1435                        // remove the ray from the leaf
1436                        // find the ray in the leaf and swap it with the last ray...
[462]1437                        VspKdLeaf *leaf = (VspKdLeaf *)data.mNode;
[465]1438
1439                        if (!leaf->Mailed())
[410]1440                        {
1441                                leaf->Mail();
[465]1442
[410]1443                                if (affectedLeaves)
1444                                        affectedLeaves->push_back(leaf);
[465]1445
1446                                if (removeAllScheduledRays)
[410]1447                                {
[420]1448                                        int tail = (int)leaf->mRays.size() - 1;
[406]1449
[465]1450                                        for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
[410]1451                                        {
[465]1452                                                if (leaf->mRays[i].mRay->ScheduledForRemoval())
[410]1453                                                {
1454                                                        // find a ray to replace it with
[465]1455                                                        while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
[410]1456                                                        {
[414]1457                                                                ++ mStat.removedRayRefs;
[420]1458                                                                leaf->mRays[tail].mRay->Unref();
1459                                                                leaf->mRays.pop_back();
[465]1460
[410]1461                                                                -- tail;
1462                                                        }
[465]1463
[410]1464                                                        if (tail < i)
1465                                                                break;
[465]1466
[414]1467                                                        ++ mStat.removedRayRefs;
[465]1468
[420]1469                                                        leaf->mRays[i].mRay->Unref();
1470                                                        leaf->mRays[i] = leaf->mRays[tail];
1471                                                        leaf->mRays.pop_back();
[410]1472                                                        -- tail;
1473                                                }
1474                                        }
1475                                }
1476                        }
[465]1477
[410]1478                        if (!removeAllScheduledRays)
[465]1479                                for (int i=0; i < (int)leaf->mRays.size(); i++)
[410]1480                                {
[465]1481                                        if (leaf->mRays[i].mRay == ray)
[410]1482                                        {
[414]1483                                                ++ mStat.removedRayRefs;
[410]1484                                                ray->Unref();
[420]1485                                                leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
1486                                                leaf->mRays.pop_back();
[410]1487                                            // check this ray again
1488                                                break;
1489                                        }
1490                                }
1491                }
[408]1492        }
1493
[465]1494        if (ray->RefCount() != 0)
[410]1495        {
1496                cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1497                exit(1);
1498        }
1499
[406]1500}
1501
[410]1502void VspKdTree::AddRay(VssRay *ray)
[406]1503{
[410]1504        stack<RayTraversalData> tstack;
[465]1505
[420]1506        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[465]1507
[410]1508        RayTraversalData data;
[465]1509
1510        while (!tstack.empty())
[410]1511        {
1512                data = tstack.top();
1513                tstack.pop();
[406]1514
[465]1515                if (!data.mNode->IsLeaf())
[410]1516                {
1517                        TraverseInternalNode(data, tstack);
[465]1518                }
1519                else
[410]1520                {
1521                        // remove the ray from the leaf
[420]1522                        // find the ray in the leaf and swap it with the last ray
[462]1523                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(data.mNode);
[423]1524
[419]1525                        leaf->AddRay(data.mRayData);
[413]1526                        ++ mStat.addedRayRefs;
[410]1527                }
1528        }
[406]1529}
1530
[465]1531void VspKdTree::TraverseInternalNode(RayTraversalData &data,
[410]1532                                                                         stack<RayTraversalData> &tstack)
[406]1533{
[462]1534        VspKdInterior *in = dynamic_cast<VspKdInterior *>(data.mNode);
[465]1535
1536        if (in->mAxis <= VspKdNode::SPLIT_Z)
[410]1537        {
1538            // determine the side of this ray with respect to the plane
[485]1539                float t;
1540                const int side = in->ComputeRayIntersection(data.mRayData, t);
[465]1541
1542                if (side == 0)
[410]1543                {
[465]1544                        if (data.mRayData.mRay->HasPosDir(in->mAxis))
[433]1545                        {
[465]1546                                tstack.push(RayTraversalData(in->GetBack(),
[463]1547                                                        RayInfo(data.mRayData.mRay, data.mRayData.mMinT, t)));
[465]1548
[433]1549                                tstack.push(RayTraversalData(in->GetFront(),
[463]1550                                                        RayInfo(data.mRayData.mRay, t, data.mRayData.mMaxT)));
[465]1551
1552                        }
1553                        else
[433]1554                        {
1555                                tstack.push(RayTraversalData(in->GetBack(),
1556                                                        RayInfo(data.mRayData.mRay,
[463]1557                                                        t,
[433]1558                                                        data.mRayData.mMaxT)));
1559                                tstack.push(RayTraversalData(in->GetFront(),
1560                                                        RayInfo(data.mRayData.mRay,
1561                                                        data.mRayData.mMinT,
[463]1562                                                        t)));
[433]1563                        }
[410]1564                }
[433]1565                else
1566                        if (side == 1)
1567                                tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1568                        else
1569                                tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[465]1570        }
1571        else
[410]1572        {
1573                // directional split
[419]1574                if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
1575                        tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
[408]1576                else
[419]1577                        tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[410]1578        }
[406]1579}
1580
[408]1581
[462]1582int VspKdTree::CollapseSubtree(VspKdNode *sroot, const int time)
[406]1583{
[410]1584        // first count all rays in the subtree
1585        // use mail 1 for this purpose
[462]1586        stack<VspKdNode *> tstack;
[465]1587
[410]1588        int rayCount = 0;
1589        int totalRayCount = 0;
1590        int collapsedNodes = 0;
[408]1591
1592#if DEBUG_COLLAPSE
[420]1593        cout << "Collapsing subtree" << endl;
[428]1594        cout << "accessTime=" << sroot->GetAccessTime() << endl;
[420]1595        cout << "depth=" << (int)sroot->depth << endl;
[408]1596#endif
1597
[410]1598        // tstat.collapsedSubtrees++;
1599        // tstat.collapseDepths += (int)sroot->depth;
1600        // tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1601        tstack.push(sroot);
1602        VssRay::NewMail();
[465]1603
1604        while (!tstack.empty())
[410]1605        {
1606                ++ collapsedNodes;
[465]1607
[462]1608                VspKdNode *node = tstack.top();
[410]1609                tstack.pop();
[465]1610                if (node->IsLeaf())
[410]1611                {
[462]1612                        VspKdLeaf *leaf = (VspKdLeaf *) node;
[465]1613
[420]1614                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1615                                        ri != leaf->mRays.end(); ++ ri)
[410]1616                        {
1617                                ++ totalRayCount;
[465]1618
1619                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed())
[410]1620                                {
[408]1621                                        (*ri).mRay->Mail();
[410]1622                                        ++ rayCount;
[408]1623                                }
[410]1624                        }
[465]1625                }
1626                else
[410]1627                {
[472]1628                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetFront());
1629                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetBack());
[410]1630                }
1631        }
[465]1632
[410]1633        VssRay::NewMail();
[406]1634
[410]1635        // create a new node that will hold the rays
[462]1636        VspKdLeaf *newLeaf = new VspKdLeaf(sroot->mParent, rayCount);
[465]1637
[500]1638        VspKdInterior *parent =
1639                dynamic_cast<VspKdInterior *>(newLeaf->mParent);
1640
1641        if (parent)
[472]1642        {
[500]1643                parent->ReplaceChildLink(sroot, newLeaf);
[472]1644        }
[410]1645        tstack.push(sroot);
[406]1646
[465]1647        while (!tstack.empty())
[410]1648        {
[462]1649                VspKdNode *node = tstack.top();
[410]1650                tstack.pop();
[408]1651
[465]1652                if (node->IsLeaf())
[410]1653                {
[462]1654                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[408]1655
[420]1656                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1657                                        ri != leaf->mRays.end(); ++ ri)
[410]1658                        {
[465]1659                                // unref this ray from the old node
1660                                if ((*ri).mRay->IsActive())
[410]1661                                {
[408]1662                                        (*ri).mRay->Unref();
[465]1663                                        if (!(*ri).mRay->Mailed())
[410]1664                                        {
[408]1665                                                (*ri).mRay->Mail();
1666                                                newLeaf->AddRay(*ri);
1667                                        }
[465]1668                                }
[410]1669                                else
[408]1670                                        (*ri).mRay->Unref();
[410]1671                        }
[465]1672                }
1673                else
[410]1674                {
[465]1675                        VspKdInterior *interior =
[462]1676                                dynamic_cast<VspKdInterior *>(node);
[419]1677                        tstack.push(interior->GetBack());
1678                        tstack.push(interior->GetFront());
[410]1679                }
1680        }
[465]1681
[410]1682        // delete the node and all its children
[420]1683        DEL_PTR(sroot);
[465]1684
[462]1685        // for(VspKdNode::SRayContainer::iterator ri = newleaf->mRays.begin();
[420]1686    //      ri != newleaf->mRays.end(); ++ ri)
[410]1687        //     (*ri).ray->UnMail(2);
[408]1688
1689#if DEBUG_COLLAPSE
[423]1690        cout << "Total memory before=" << GetMemUsage() << endl;
[408]1691#endif
1692
[414]1693        mStat.nodes -= collapsedNodes - 1;
1694        mStat.rayRefs -= totalRayCount - rayCount;
[465]1695
[408]1696#if DEBUG_COLLAPSE
[410]1697        cout << "collapsed nodes" << collapsedNodes << endl;
1698        cout << "collapsed rays" << totalRayCount - rayCount << endl;
1699        cout << "Total memory after=" << GetMemUsage() << endl;
1700        cout << "================================" << endl;
[408]1701#endif
1702
1703        //  tstat.collapsedNodes += collapsedNodes;
1704        //  tstat.collapsedRays += totalRayCount - rayCount;
[410]1705    return totalRayCount - rayCount;
[406]1706}
1707
[408]1708
[462]1709int VspKdTree::GetPvsSize(VspKdNode *node, const AxisAlignedBox3 &box) const
[406]1710{
[462]1711        stack<VspKdNode *> tstack;
[419]1712        tstack.push(mRoot);
[408]1713
1714        Intersectable::NewMail();
1715        int pvsSize = 0;
[465]1716
1717        while (!tstack.empty())
[410]1718        {
[462]1719                VspKdNode *node = tstack.top();
[410]1720                tstack.pop();
[465]1721
1722          if (node->IsLeaf())
[410]1723                {
[462]1724                        VspKdLeaf *leaf = (VspKdLeaf *)node;
[420]1725                        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1726                                 ri != leaf->mRays.end(); ++ ri)
[410]1727                        {
[465]1728                                if ((*ri).mRay->IsActive())
[410]1729                                {
[408]1730                                        Intersectable *object;
1731#if BIDIRECTIONAL_RAY
1732                                        object = (*ri).mRay->mOriginObject;
[465]1733
1734                                        if (object && !object->Mailed())
[410]1735                                        {
1736                                                ++ pvsSize;
[408]1737                                                object->Mail();
1738                                        }
1739#endif
1740                                        object = (*ri).mRay->mTerminationObject;
[465]1741                                        if (object && !object->Mailed())
[410]1742                                        {
1743                                                ++ pvsSize;
[408]1744                                                object->Mail();
1745                                        }
1746                                }
[465]1747                        }
[410]1748                }
[465]1749                else
[410]1750                {
[462]1751                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
[465]1752
1753                        if (in->mAxis < 3)
[410]1754                        {
[419]1755                                if (box.Max(in->mAxis) >= in->mPosition)
1756                                        tstack.push(in->GetFront());
[465]1757
[419]1758                                if (box.Min(in->mAxis) <= in->mPosition)
1759                                        tstack.push(in->GetBack());
[465]1760                        }
1761                        else
[410]1762                        {
[408]1763                                // both nodes for directional splits
[419]1764                                tstack.push(in->GetFront());
1765                                tstack.push(in->GetBack());
[408]1766                        }
1767                }
1768        }
[410]1769
[408]1770        return pvsSize;
[406]1771}
1772
[410]1773void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1774                                                                                         float &maxRayContribution,
1775                                                                                         float &avgRayContribution)
[406]1776{
[462]1777        stack<VspKdNode *> tstack;
[419]1778        tstack.push(mRoot);
[465]1779
[408]1780        minRayContribution = 1.0f;
1781        maxRayContribution = 0.0f;
[410]1782
[408]1783        float sumRayContribution = 0.0f;
1784        int leaves = 0;
[406]1785
[465]1786        while (!tstack.empty())
[410]1787        {
[462]1788                VspKdNode *node = tstack.top();
[410]1789                tstack.pop();
[465]1790                if (node->IsLeaf())
[410]1791                {
[408]1792                        leaves++;
[462]1793                        VspKdLeaf *leaf = (VspKdLeaf *)node;
[408]1794                        float c = leaf->GetAvgRayContribution();
[410]1795
[408]1796                        if (c > maxRayContribution)
1797                                maxRayContribution = c;
1798                        if (c < minRayContribution)
1799                                minRayContribution = c;
1800                        sumRayContribution += c;
[465]1801
1802                }
[482]1803                else
1804                {
[462]1805                        VspKdInterior *in = (VspKdInterior *)node;
[482]1806                       
[419]1807                        tstack.push(in->GetFront());
1808                        tstack.push(in->GetBack());
[408]1809                }
1810        }
[465]1811
[423]1812        Debug << "sum=" << sumRayContribution << endl;
1813        Debug << "leaves=" << leaves << endl;
1814        avgRayContribution = sumRayContribution / (float)leaves;
[406]1815}
1816
1817
[410]1818int VspKdTree::GenerateRays(const float ratioPerLeaf,
1819                                                        SimpleRayContainer &rays)
[406]1820{
[462]1821        stack<VspKdNode *> tstack;
[419]1822        tstack.push(mRoot);
[465]1823
1824        while (!tstack.empty())
[410]1825        {
[462]1826                VspKdNode *node = tstack.top();
[410]1827                tstack.pop();
[465]1828
1829                if (node->IsLeaf())
[410]1830                {
[462]1831                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[423]1832
[408]1833                        float c = leaf->GetAvgRayContribution();
[411]1834                        int num = (int)(c*ratioPerLeaf + 0.5);
[408]1835                        //                      cout<<num<<" ";
1836
[465]1837                        for (int i=0; i < num; i++)
[410]1838                        {
[408]1839                                Vector3 origin = GetBBox(leaf).GetRandomPoint();
[419]1840                                /*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
[408]1841                                Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1842                                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[419]1843                                rays.push_back(SimpleRay(origin, direction));*/
[408]1844                        }
[465]1845
1846                }
1847                else
[410]1848                {
[465]1849                        VspKdInterior *in =
[462]1850                                dynamic_cast<VspKdInterior *>(node);
[482]1851               
[419]1852                        tstack.push(in->GetFront());
1853                        tstack.push(in->GetBack());
[408]1854                }
[406]1855        }
1856
[411]1857        return (int)rays.size();
[406]1858}
1859
1860
[410]1861float VspKdTree::GetAvgPvsSize()
[406]1862{
[462]1863        stack<VspKdNode *> tstack;
[419]1864        tstack.push(mRoot);
[406]1865
[408]1866        int sumPvs = 0;
1867        int leaves = 0;
[465]1868
1869        while (!tstack.empty())
[410]1870        {
[462]1871                VspKdNode *node = tstack.top();
[410]1872                tstack.pop();
[465]1873
1874                if (node->IsLeaf())
[410]1875                {
[462]1876                        VspKdLeaf *leaf = (VspKdLeaf *)node;
[408]1877                        // update pvs size
1878                        leaf->UpdatePvsSize();
1879                        sumPvs += leaf->GetPvsSize();
1880                        leaves++;
[465]1881                }
1882                else
[410]1883                {
[462]1884                        VspKdInterior *in = (VspKdInterior *)node;
[482]1885                       
[419]1886                        tstack.push(in->GetFront());
1887                        tstack.push(in->GetBack());
[406]1888                }
1889        }
[408]1890
[410]1891        return sumPvs / (float)leaves;
[418]1892}
1893
[462]1894VspKdNode *VspKdTree::GetRoot() const
[418]1895{
1896        return mRoot;
1897}
1898
[462]1899AxisAlignedBox3 VspKdTree::GetBBox(VspKdNode *node) const
[418]1900{
[419]1901        if (node->mParent == NULL)
1902                return mBox;
[418]1903
1904        if (!node->IsLeaf())
[500]1905        {
1906                if (node->Type() == VspKdNode::EIntermediate)
1907            return (dynamic_cast<VspKdInterior *>(node))->mBox;
1908                else
1909                        return (dynamic_cast<VspKdIntermediate *>(node))->mBox;
1910        }
[418]1911
[500]1912        VspKdInterior *parent = dynamic_cast<VspKdInterior *>(node->mParent);
[465]1913
[500]1914        if (parent->mAxis >= 3)
1915                return parent->mBox;
1916
1917        AxisAlignedBox3 box(parent->mBox);
1918        if (parent->GetFront() == node)
1919                box.SetMin(parent->mAxis, parent->mPosition);
[418]1920        else
[500]1921                box.SetMax(parent->mAxis, parent->mPosition);
[419]1922
1923        return box;
[418]1924}
1925
[465]1926int     VspKdTree::GetRootPvsSize() const
[418]1927{
1928        return GetPvsSize(mRoot, mBox);
1929}
[419]1930
[465]1931const VspKdStatistics &VspKdTree::GetStatistics() const
[419]1932{
1933        return mStat;
[420]1934}
1935
1936void VspKdTree::AddRays(VssRayContainer &add)
1937{
1938        VssRayContainer remove;
1939        UpdateRays(remove, add);
1940}
[465]1941
[420]1942// return memory usage in MB
[465]1943float VspKdTree::GetMemUsage() const
[420]1944{
[465]1945        return
1946                (sizeof(VspKdTree) +
[462]1947                 mStat.Leaves() * sizeof(VspKdLeaf) +
1948                 mStat.Interior() * sizeof(VspKdInterior) +
[420]1949                 mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
1950}
[465]1951
1952float VspKdTree::GetRayMemUsage() const
[420]1953{
1954        return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
[425]1955}
1956
[482]1957
[465]1958int VspKdTree::ComputePvsSize(VspKdNode *node,
[426]1959                                                          const RayInfoContainer &globalRays) const
[425]1960{
1961        int pvsSize = 0;
1962
1963        const AxisAlignedBox3 box = GetBBox(node);
1964
1965        RayInfoContainer::const_iterator it, it_end = globalRays.end();
1966
[495]1967        Intersectable::NewMail();
[482]1968
[426]1969        // warning: implicit conversion from VssRay to Ray
[425]1970        for (it = globalRays.begin(); it != globalRays.end(); ++ it)
[495]1971        {
1972                VssRay *ray = (*it).mRay;
[465]1973
[495]1974                // ray intersects node bounding box
1975                if (box.GetMinMaxT(*ray, NULL, NULL))
1976                {
1977                        if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
1978                        {
1979                                ray->mTerminationObject->Mail();
1980                                ++ pvsSize;
1981                        }
1982                }
1983        }
[425]1984        return pvsSize;
[428]1985}
1986
[482]1987
[462]1988void VspKdTree::CollectLeaves(vector<VspKdLeaf *> &leaves) const
[428]1989{
[462]1990        stack<VspKdNode *> nodeStack;
[428]1991        nodeStack.push(mRoot);
1992
[465]1993        while (!nodeStack.empty())
[428]1994        {
[462]1995                VspKdNode *node = nodeStack.top();
[465]1996
[428]1997                nodeStack.pop();
[465]1998
1999                if (node->IsLeaf())
[428]2000                {
[462]2001                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[428]2002                        leaves.push_back(leaf);
[465]2003                }
2004                else
[428]2005                {
[462]2006                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
[428]2007
[448]2008                        nodeStack.push(interior->GetBack());
2009                        nodeStack.push(interior->GetFront());
[428]2010                }
2011        }
[452]2012}
2013
[503]2014
[462]2015int VspKdTree::FindNeighbors(VspKdLeaf *n,
2016                                                         vector<VspKdLeaf *> &neighbors,
[452]2017                                                         bool onlyUnmailed)
2018{
[462]2019        stack<VspKdNode *> nodeStack;
[465]2020
[452]2021        nodeStack.push(mRoot);
2022
2023        AxisAlignedBox3 box = GetBBox(n);
2024
[465]2025        while (!nodeStack.empty())
[452]2026        {
[462]2027                VspKdNode *node = nodeStack.top();
[452]2028                nodeStack.pop();
2029
[465]2030                if (node->IsLeaf())
[452]2031                {
[462]2032                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[452]2033
[465]2034                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
[452]2035                                neighbors.push_back(leaf);
[465]2036                }
2037                else
[452]2038                {
[462]2039                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
[452]2040
2041                        if (interior->mPosition > box.Max(interior->mAxis))
2042                                nodeStack.push(interior->mBack);
2043                        else
2044                        {
2045                                if (interior->mPosition < box.Min(interior->mAxis))
2046                                        nodeStack.push(interior->mFront);
[465]2047                                else
[452]2048                                {
2049                                        // random decision
2050                                        nodeStack.push(interior->mBack);
2051                                        nodeStack.push(interior->mFront);
2052                                }
2053                        }
2054                }
2055        }
2056
2057        return (int)neighbors.size();
2058}
2059
[462]2060
2061void VspKdTree::SetViewCellsManager(ViewCellsManager *vcm)
[452]2062{
[462]2063        mViewCellsManager = vcm;
2064}
[452]2065
[462]2066
[468]2067int VspKdTree::CastLineSegment(const Vector3 &origin,
2068                                                           const Vector3 &termination,
2069                                                           vector<ViewCell *> &viewcells)
[462]2070{
2071        int hits = 0;
[465]2072
[468]2073        float mint = 0.0f, maxt = 1.0f;
2074        const Vector3 dir = termination - origin;
[465]2075
[483]2076        stack<LineTraversalData> tStack;
[462]2077
2078        Intersectable::NewMail();
2079
[468]2080        Vector3 entp = origin;
2081        Vector3 extp = termination;
[465]2082
[462]2083        VspKdNode *node = mRoot;
2084        VspKdNode *farChild;
2085
2086        float position;
2087        int axis;
[465]2088
2089        while (1)
[462]2090        {
[465]2091                if (!node->IsLeaf())
[462]2092                {
2093                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
2094                        position = in->mPosition;
2095                        axis = in->mAxis;
2096
[465]2097                        if (entp[axis] <= position)
[462]2098                        {
[465]2099                                if (extp[axis] <= position)
[462]2100                                {
2101                                        node = in->mBack;
2102                                        // cases N1,N2,N3,P5,Z2,Z3
2103                                        continue;
[465]2104                                } else
[462]2105                                {
2106                                        // case N4
2107                                        node = in->mBack;
2108                                        farChild = in->mFront;
2109                                }
2110                        }
[465]2111                        else
[462]2112                        {
[465]2113                                if (position <= extp[axis])
[462]2114                                {
2115                                        node = in->mFront;
2116                                        // cases P1,P2,P3,N5,Z1
2117                                        continue;
[465]2118                                }
2119                                else
[462]2120                                {
2121                                        node = in->mFront;
2122                                        farChild = in->mBack;
2123                                        // case P4
2124                                }
2125                        }
2126
2127                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2128                        // case N4 or P4
[468]2129                        float tdist = (position - origin[axis]) / dir[axis];
[483]2130                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
[468]2131                        extp = origin + dir * tdist;
[462]2132                        maxt = tdist;
[465]2133                }
2134                else
[462]2135                {
2136                        // compute intersection with all objects in this leaf
[468]2137                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2138                        ViewCell *vc = leaf->GetViewCell();
[465]2139
[468]2140                        if (!vc->Mailed())
[462]2141                        {
[468]2142                                vc->Mail();
2143                                viewcells.push_back(vc);
[469]2144                                ++ hits;
[462]2145                        }
[465]2146
[485]2147                        if (0) //matt TODO: REMOVE LATER
[483]2148                                leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2149
[462]2150                        // get the next node from the stack
2151                        if (tStack.empty())
2152                                break;
[465]2153
[462]2154                        entp = extp;
2155                        mint = maxt;
[468]2156                       
[483]2157                        LineTraversalData &s  = tStack.top();
[462]2158                        node = s.mNode;
2159                        extp = s.mExitPoint;
2160                        maxt = s.mMaxT;
2161                        tStack.pop();
2162                }
2163        }
[465]2164
[462]2165        return hits;
2166}
2167
2168
2169void VspKdTree::GenerateViewCell(VspKdLeaf *leaf)
2170{
2171        RayInfoContainer::const_iterator it, it_end = leaf->GetRays().end();
2172
2173        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(mViewCellsManager->GenerateViewCell());
2174        leaf->SetViewCell(vc);
2175
[479]2176        vc->SetVolume(GetBBox(leaf).GetVolume());
2177        vc->SetArea(GetBBox(leaf).SurfaceArea());
2178
[580]2179        vc->mLeaf = leaf;
[462]2180
2181        for (it = leaf->GetRays().begin(); it != it_end; ++ it)
2182        {
2183                VssRay *ray = (*it).mRay;
2184
[463]2185                if (ray->mTerminationObject)
[556]2186                  vc->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf);
[465]2187
[463]2188                if (ray->mOriginObject)
[556]2189                  vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf);
[462]2190        }
2191}
2192
2193
[589]2194bool VspKdTree::MergeViewCells(VspKdLeaf *l1, VspKdLeaf *l2)
2195{
2196        //-- change pointer to view cells of all leaves associated
2197        //-- with the previous view cells
2198        VspKdViewCell *fVc = l1->GetViewCell();
2199        VspKdViewCell *bVc = l2->GetViewCell();
2200
2201        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(
2202                mViewCellsManager->MergeViewCells(fVc, bVc));
2203
2204        // if merge was unsuccessful
2205        if (!vc) return false;
2206// TODO
[728]2207#if HAS_TO_BE_REDONE
[589]2208        // set new size of view cell
[479]2209        vc->SetVolume(fVc->GetVolume() + bVc->GetVolume());
[482]2210        // new area
[589]2211        vc->SetArea(fVc->GetArea() + bVc->GetArea());
2212
2213        vector<VspKdLeaf *> fLeaves = fVc->mLeaves;
2214        vector<VspKdLeaf *> bLeaves = bVc->mLeaves;
2215
2216        vector<VspKdLeaf *>::const_iterator it;
2217
2218        //-- change view cells of all the other leaves the view cell belongs to
2219        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
2220        {
2221                (*it)->SetViewCell(vc);
2222                vc->mLeaves.push_back(*it);
2223        }
2224
2225        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
2226        {
2227                (*it)->SetViewCell(vc);
2228                vc->mLeaves.push_back(*it);
2229        }
2230#endif
[465]2231        vc->Mail();
[589]2232
2233        // clean up old view cells
2234        DEL_PTR(fVc);
2235        DEL_PTR(bVc);
2236
[463]2237        return true;
[462]2238}
2239
[483]2240void VspKdTree::CollectMergeCandidates()
[462]2241{
[580]2242        VspKdMergeCandidate::sOverallCost = 0;
[462]2243        vector<VspKdLeaf *> leaves;
[465]2244
[462]2245        // collect the leaves, e.g., the "voxels" that will build the view cells
[452]2246        CollectLeaves(leaves);
2247
[462]2248        VspKdLeaf::NewMail();
[452]2249
[462]2250        vector<VspKdLeaf *>::const_iterator it, it_end = leaves.end();
[452]2251
[462]2252        // generate view cells
[452]2253        for (it = leaves.begin(); it != it_end; ++ it)
[462]2254                GenerateViewCell(*it);
2255
2256        // find merge candidates and push them into queue
2257        for (it = leaves.begin(); it != it_end; ++ it)
[452]2258        {
[462]2259                VspKdLeaf *leaf = *it;
[485]2260                ViewCell *vc = leaf->GetViewCell();
[462]2261                // no leaf is part of two merge candidates
[485]2262                leaf->Mail();
[580]2263                VspKdMergeCandidate::sOverallCost +=
[485]2264                        vc->GetVolume() * vc->GetPvs().GetSize();
2265                vector<VspKdLeaf *> neighbors;
2266                FindNeighbors(leaf, neighbors, true);
[462]2267
[485]2268                vector<VspKdLeaf *>::const_iterator nit,
2269                        nit_end = neighbors.end();
[452]2270
[485]2271                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
2272                {
[580]2273                        mMergeQueue.push(VspKdMergeCandidate(leaf, *nit));
[462]2274                }
2275        }
[483]2276}
2277
2278
[485]2279int VspKdTree::MergeViewCells(const VssRayContainer &rays)
[483]2280{
[485]2281        CollectMergeCandidates();
2282
[462]2283        int merged = 0;
2284
[485]2285        int vcSize = mStat.Leaves();
[462]2286        // use priority queue to merge leaves
[483]2287        while (!mMergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
2288                   (mMergeQueue.top().GetMergeCost() <
[580]2289                    mMergeMaxCostRatio * VspKdMergeCandidate::sOverallCost))
[462]2290        {
[485]2291                //Debug << "mergecost: " << mergeQueue.top().GetMergeCost() /
[580]2292                //VspKdMergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
2293                VspKdMergeCandidate mc = mMergeQueue.top();
[483]2294                mMergeQueue.pop();
[465]2295
2296                // both view cells equal!
[462]2297                if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
2298                        continue;
2299
2300                if (mc.Valid())
2301                {
[465]2302                        ViewCell::NewMail();
2303                        MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
[462]2304
2305                        ++ merged;
2306                        -- vcSize;
[472]2307                        // increase absolute merge cost
[580]2308                        VspKdMergeCandidate::sOverallCost += mc.GetMergeCost();
[452]2309                }
[462]2310                // merge candidate not valid, because one of the leaves was already
2311                // merged with another one
[465]2312                else
[462]2313                {
2314                        // validate and reinsert into queue
2315                        mc.SetValid();
[483]2316                        mMergeQueue.push(mc);
[465]2317                        //Debug << "validate " << mc.GetMergeCost() << endl;
[462]2318                }
[452]2319        }
[462]2320
[485]2321        Debug << "merged " << merged << " of " << mStat.Leaves() << " leaves" << endl;
[465]2322
[469]2323        //TODO: should return sample contributions
[462]2324        return merged;
[452]2325}
2326
2327
[517]2328void VspKdTree::RepairViewCellsLeafLists()
[589]2329{
[465]2330        // list not valid anymore => clear
[589]2331        stack<VspKdNode *> nodeStack;
2332        nodeStack.push(mRoot);
2333
2334        ViewCell::NewMail();
2335
2336        while (!nodeStack.empty())
2337        {
2338                VspKdNode *node = nodeStack.top();
2339                nodeStack.pop();
2340
2341                if (node->IsLeaf())
2342                {
2343                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2344
2345                        VspKdViewCell *viewCell = leaf->GetViewCell();
[728]2346#if HAS_TO_BE_REDONE
[589]2347                        if (!viewCell->Mailed())
2348                        {
2349                                viewCell->mLeaves.clear();
2350                                viewCell->Mail();
2351                        }
2352
2353                        viewCell->mLeaves.push_back(leaf);
2354#endif
2355                }
2356                else
2357                {
2358                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2359
2360                        nodeStack.push(interior->GetFront());
2361                        nodeStack.push(interior->GetBack());
2362                }
[465]2363        }
2364}
2365
[479]2366
[501]2367VspKdNode *VspKdTree::CollapseTree(VspKdNode *node, int &collapsed)
[465]2368{
[501]2369        if (node->IsLeaf())
[465]2370                return node;
2371
[501]2372        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2373
2374        VspKdNode *front = CollapseTree(interior->GetFront(), collapsed);
2375        VspKdNode *back = CollapseTree(interior->GetBack(), collapsed);
2376
[465]2377        if (front->IsLeaf() && back->IsLeaf())
2378        {
2379                VspKdLeaf *frontLeaf = dynamic_cast<VspKdLeaf *>(front);
2380                VspKdLeaf *backLeaf = dynamic_cast<VspKdLeaf *>(back);
2381
2382                //-- collapse tree
2383                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
2384                {
2385                        VspKdViewCell *vc = frontLeaf->GetViewCell();
2386
2387                        VspKdLeaf *leaf = new VspKdLeaf(interior->GetParent(), 0);
2388                        leaf->SetViewCell(vc);
[501]2389                        //leaf->SetTreeValid(frontLeaf->TreeValid());
[465]2390
2391                        // replace a link from node's parent
[501]2392                        if (leaf->mParent)
2393                        {
2394                                // replace a link from node's parent
2395                                VspKdInterior *parent =
2396                                        dynamic_cast<VspKdInterior *>(leaf->mParent);
[500]2397                                parent->ReplaceChildLink(node, leaf);
[501]2398                        }
2399                        ++ collapsed;
[465]2400                        delete interior;
2401
2402                        return leaf;
2403                }
2404        }
2405
[495]2406        return node;
2407}
2408
2409
[501]2410int VspKdTree::CollapseTree()
[495]2411{
[501]2412        int collapsed = 0;
2413        CollapseTree(mRoot, collapsed);
[485]2414        // revalidate leaves
[517]2415        RepairViewCellsLeafLists();
[501]2416
2417        return collapsed;
[465]2418}
2419
[501]2420
[486]2421int VspKdTree::RefineViewCells(const VssRayContainer &rays)
[473]2422{
[486]2423        int shuffled = 0;
2424
2425        Debug << "refining " << (int)mMergeQueue.size() << " candidates " << endl;
2426        BspLeaf::NewMail();
2427
2428        // Use priority queue of remaining leaf pairs
2429        // These candidates either share the same view cells or
2430        // are border leaves which share a boundary.
2431        // We test if they can be shuffled, i.e.,
2432        // either one leaf is made part of one view cell or the other
2433        // leaf is made part of the other view cell. It is tested if the
2434        // remaining view cells are "better" than the old ones.
2435        while (!mMergeQueue.empty())
2436        {
[580]2437                VspKdMergeCandidate mc = mMergeQueue.top();
[486]2438                mMergeQueue.pop();
2439
2440                // both view cells equal or already shuffled
2441                if ((mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell()) ||
2442                        (mc.GetLeaf1()->Mailed()) || (mc.GetLeaf2()->Mailed()))
2443                        continue;
2444               
2445                // candidate for shuffling
2446                const bool wasShuffled = false;
2447                        //ShuffleLeaves(mc.GetLeaf1(), mc.GetLeaf2());
[495]2448                        // matt: restore
[486]2449                //-- stats
2450                if (wasShuffled)
2451                        ++ shuffled;
2452        }
2453
2454        return shuffled;
[473]2455}
2456
2457
[486]2458inline int AddedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2459{
2460        return pvs1.AddPvs(pvs2);
2461}
2462
2463
2464inline int SubtractedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2465{
2466        return pvs1.SubtractPvs(pvs2);
2467}
2468
2469
[580]2470float EvalShuffleCost(VspKdLeaf *leaf, VspKdViewCell *vc1, VspKdViewCell *vc2)
[486]2471{
2472        const int pvs1 = SubtractedPvsSize(vc1->GetPvs(), *leaf->mPvs);
2473        const int pvs2 = AddedPvsSize(vc2->GetPvs(), *leaf->mPvs);
2474
[495]2475        const float vol1 = vc1->GetVolume() - leaf->mVolume;
2476        const float vol2 = vc2->GetVolume() + leaf->mVolume;
[486]2477
[495]2478        const float cost1 = pvs1 * vol1;
2479        const float cost2 = pvs2 * vol2;
[486]2480
2481        return cost1 + cost2;
2482}
2483
2484
2485void VspKdTree::ShuffleLeaf(VspKdLeaf *leaf,
2486                                                        VspKdViewCell *vc1,
2487                                                        VspKdViewCell *vc2) const
2488{
2489        // compute new pvs and area
[728]2490#if HAS_TO_BE_REDONE
[486]2491        vc1->GetPvs().SubtractPvs(*leaf->mPvs);
2492        vc2->GetPvs().AddPvs(*leaf->mPvs);
2493       
[495]2494        vc1->SetArea(vc1->GetVolume() - leaf->mVolume);
2495        vc2->SetArea(vc2->GetVolume() + leaf->mVolume);
[486]2496
2497        /// add to second view cell
2498        vc2->mLeaves.push_back(leaf);
2499
2500        // erase leaf from old view cell
2501        vector<VspKdLeaf *>::iterator it = vc1->mLeaves.begin();
2502
2503        for (; *it != leaf; ++ it);
2504        vc1->mLeaves.erase(it);
2505
2506        leaf->SetViewCell(vc2);  // finally change view cell
[580]2507#endif
[486]2508}
2509
2510
2511bool VspKdTree::ShuffleLeaves(VspKdLeaf *leaf1, VspKdLeaf *leaf2) const
2512{
2513        VspKdViewCell *vc1 = leaf1->GetViewCell();
2514        VspKdViewCell *vc2 = leaf2->GetViewCell();
2515
2516        const float cost1 = vc1->GetPvs().GetSize() * vc1->GetArea();
2517        const float cost2 = vc2->GetPvs().GetSize() * vc2->GetArea();
2518
2519        const float oldCost = cost1 + cost2;
2520       
2521        float shuffledCost1 = Limits::Infinity;
2522        float shuffledCost2 = Limits::Infinity;
2523
2524        // the view cell should not be empty after the shuffle
[580]2525//todo
2526        /*      if (vc1->mLeaves.size() > 1)
2527                shuffledCost1 = EvalShuffleCost(leaf1, vc1, vc2);
[486]2528        if (vc2->mLeaves.size() > 1)
[580]2529                shuffledCost2 = EvalShuffleCost(leaf2, vc2, vc1);
[486]2530
[580]2531*/      // shuffling unsuccessful
[486]2532        if ((oldCost <= shuffledCost1) && (oldCost <= shuffledCost2))
2533                return false;
2534       
2535        if (shuffledCost1 < shuffledCost2)
2536        {
2537                ShuffleLeaf(leaf1, vc1, vc2);
2538                leaf1->Mail();
2539        }
2540        else
2541        {
2542                ShuffleLeaf(leaf2, vc2, vc1);
2543                leaf2->Mail();
2544        }
2545
2546        return true;
2547}
2548
2549
2550
2551
[452]2552/*********************************************************************/
[580]2553/*                VspKdMergeCandidate implementation                      */
[452]2554/*********************************************************************/
2555
2556
[580]2557VspKdMergeCandidate::VspKdMergeCandidate(VspKdLeaf *l1, VspKdLeaf *l2):
[462]2558mMergeCost(0),
2559mLeaf1(l1),
2560mLeaf2(l2),
2561mLeaf1Id(l1->GetViewCell()->mMailbox),
2562mLeaf2Id(l2->GetViewCell()->mMailbox)
2563{
2564        EvalMergeCost();
2565}
[452]2566
[580]2567float VspKdMergeCandidate::GetLeaf1Cost() const
[472]2568{
2569        VspKdViewCell *vc = mLeaf1->GetViewCell();
2570        return vc->GetPvs().GetSize() * vc->GetVolume();
2571}
[462]2572
[580]2573float VspKdMergeCandidate::GetLeaf2Cost() const
[472]2574{
2575        VspKdViewCell *vc = mLeaf2->GetViewCell();
2576        return vc->GetPvs().GetSize() * vc->GetVolume();
2577}
2578
[580]2579void VspKdMergeCandidate::EvalMergeCost()
[452]2580{
[465]2581        //-- compute pvs difference
[463]2582        VspKdViewCell *vc1 = mLeaf1->GetViewCell();
[462]2583        VspKdViewCell *vc2 = mLeaf2->GetViewCell();
[452]2584
[463]2585        const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
[462]2586        const int vcPvs = diff1 + vc1->GetPvs().GetSize();
[453]2587
[465]2588        //-- compute ratio of old cost
[462]2589        //-- (added size of left and right view cell times pvs size)
2590        //-- to new rendering cost (size of merged view cell times pvs size)
[472]2591        const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();
2592       
[465]2593        const float newCost =
[469]2594                (float)vcPvs * (vc1->GetVolume() + vc2->GetVolume());
[453]2595
[465]2596        mMergeCost = newCost - oldCost;
[462]2597
[465]2598//      if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
2599//              mMergeCost += 1.0;
[462]2600}
2601
[580]2602void VspKdMergeCandidate::SetLeaf1(VspKdLeaf *l)
[462]2603{
2604        mLeaf1 = l;
2605}
2606
[580]2607void VspKdMergeCandidate::SetLeaf2(VspKdLeaf *l)
[462]2608{
2609        mLeaf2 = l;
2610}
2611
[482]2612
[580]2613VspKdLeaf *VspKdMergeCandidate::GetLeaf1()
[462]2614{
2615        return mLeaf1;
2616}
2617
[482]2618
[580]2619VspKdLeaf *VspKdMergeCandidate::GetLeaf2()
[462]2620{
2621        return mLeaf2;
2622}
2623
[482]2624
[580]2625bool VspKdMergeCandidate::Valid() const
[465]2626{
2627        return
2628                (mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
[462]2629                (mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
2630}
2631
[482]2632
[580]2633float VspKdMergeCandidate::GetMergeCost() const
[462]2634{
2635        return mMergeCost;
2636}
2637
[482]2638
[580]2639void VspKdMergeCandidate::SetValid()
[462]2640{
2641        mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
2642        mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;
[465]2643
[462]2644        EvalMergeCost();
[463]2645}
Note: See TracBrowser for help on using the repository browser.