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

Revision 667, 59.5 KB checked in by mattausch, 19 years ago (diff)

test version
build script for testing. code is set uo for testing also

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)
[413]1127                ++ mStat.maxDepthNodes;
[465]1128
[421]1129        if (leaf->GetPvsSize() < mTermMinPvs)
[413]1130                ++ mStat.minPvsNodes;
[408]1131
[482]1132        mStat.accPvsSize += leaf->GetPvsSize();
1133
[426]1134        if ((int)leaf->GetRays().size() < mTermMinRays)
[413]1135                ++ mStat.minRaysNodes;
[408]1136
[426]1137        if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
[413]1138                ++ mStat.maxRayContribNodes;
[465]1139
1140        if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
[413]1141                ++ mStat.minSizeNodes;
[408]1142}
1143
1144
[465]1145inline bool VspKdTree::TerminationCriteriaMet(const VspKdLeaf *leaf,
[421]1146                                                                                          const AxisAlignedBox3 &box) const
1147{
[465]1148        return ((leaf->GetPvsSize() < mTermMinPvs) ||
[426]1149                    (leaf->mRays.size() < mTermMinRays) ||
[482]1150                        (leaf->GetAvgRayContribution() > mTermMaxRayContribution ) ||
[465]1151                        (leaf->mDepth >= mTermMaxDepth) ||
[587]1152                        (mStat.Leaves() >= mMaxViewCells) ||
[423]1153                        (SqrMagnitude(box.Size()) <= mTermMinSize));
[421]1154}
[408]1155
[421]1156
[462]1157VspKdNode *VspKdTree::SubdivideNode(VspKdLeaf *leaf,
[472]1158                                                                        const AxisAlignedBox3 &box,
1159                                                                        AxisAlignedBox3 &backBBox,
1160                                                                        AxisAlignedBox3 &frontBBox)
[408]1161{
[421]1162        if (TerminationCriteriaMet(leaf, box))
[410]1163        {
[465]1164                if (1 && leaf->mDepth >= mTermMaxDepth)
[422]1165                {
[428]1166                        Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
1167                        Debug << "Bbox: " << GetBBox(leaf) << endl;
[408]1168                }
[428]1169                //Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
[465]1170
[408]1171                return leaf;
1172        }
[465]1173
[408]1174        float position;
1175        // first count ray sides
[410]1176        int raysBack;
1177        int raysFront;
[408]1178        int pvsBack;
1179        int pvsFront;
[465]1180
[410]1181        // select subdivision axis
[425]1182        const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
[426]1183        //Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" <<       pvsFront << endl;
[408]1184
[465]1185        if (axis == -1)
[482]1186                return leaf;
1187       
[413]1188        mStat.nodes += 2;
[482]1189        ++ mStat.splits[axis];
[408]1190
[410]1191        // add the new nodes to the tree
[500]1192        VspKdInterior *node =
1193                new VspKdInterior(dynamic_cast<VspKdInterior *>(leaf->mParent));
[406]1194
[419]1195        node->mAxis = axis;
1196        node->mPosition = position;
1197        node->mBox = box;
[465]1198
[410]1199        backBBox = box;
1200        frontBBox = box;
[406]1201
[482]1202        VspKdLeaf *back = new VspKdLeaf(node, raysBack, leaf->GetMaxCostMisses());
[408]1203        back->SetPvsSize(pvsBack);
[482]1204        VspKdLeaf *front = new VspKdLeaf(node, raysFront, leaf->GetMaxCostMisses());
[408]1205        front->SetPvsSize(pvsFront);
[406]1206
[410]1207        // replace a link from node's parent
[419]1208        if (leaf->mParent)
[500]1209        {
1210                VspKdInterior *parent = dynamic_cast<VspKdInterior *>(leaf->mParent);
1211                parent->ReplaceChildLink(leaf, node);
1212        }
[501]1213
[410]1214        // and setup child links
1215        node->SetupChildLinks(back, front);
[465]1216
1217        if (axis <= VspKdNode::SPLIT_Z)
[410]1218        {
[408]1219                backBBox.SetMax(axis, position);
[410]1220                frontBBox.SetMin(axis, position);
[465]1221
[420]1222                for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1223                                ri != leaf->mRays.end(); ++ ri)
[410]1224                {
[465]1225                        if ((*ri).mRay->IsActive())
[410]1226                        {
[408]1227                                // first unref ray from the former leaf
1228                                (*ri).mRay->Unref();
[463]1229                                  float t;
[408]1230                                // determine the side of this ray with respect to the plane
[463]1231                                int side = node->ComputeRayIntersection(*ri, t);
[408]1232
[465]1233                                if (side == 0)
[410]1234                                {
[465]1235                                        if ((*ri).mRay->HasPosDir(axis))
[410]1236                                        {
[420]1237                                                back->AddRay(RayInfo((*ri).mRay,
1238                                                                                         (*ri).mMinT,
[463]1239                                                                                         t));
[420]1240                                                front->AddRay(RayInfo((*ri).mRay,
[463]1241                                                                                          t,
[420]1242                                                                                          (*ri).mMaxT));
[465]1243                                        }
1244                                        else
[410]1245                                        {
[420]1246                                                back->AddRay(RayInfo((*ri).mRay,
[463]1247                                                                                         t,
[420]1248                                                                                         (*ri).mMaxT));
1249                                                front->AddRay(RayInfo((*ri).mRay,
1250                                                                                          (*ri).mMinT,
[463]1251                                                                                          t));
[408]1252                                        }
[465]1253                                }
[410]1254                                else
[465]1255                                        if (side == 1)
[408]1256                                                front->AddRay(*ri);
1257                                        else
1258                                                back->AddRay(*ri);
[410]1259                        } else
[408]1260                                (*ri).mRay->Unref();
[410]1261                }
[465]1262        }
1263        else
[410]1264        {
1265                // rays front/back
[465]1266
[420]1267        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1268                        ri != leaf->mRays.end(); ++ ri)
[410]1269                {
[465]1270                        if ((*ri).mRay->IsActive())
[410]1271                        {
[408]1272                                // first unref ray from the former leaf
1273                                (*ri).mRay->Unref();
1274
1275                                int side;
1276                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1277                                        side = 1;
1278                                else
1279                                        side = -1;
[465]1280
[408]1281                                if (side == 1)
1282                                        front->AddRay(*ri);
1283                                else
[465]1284                                        back->AddRay(*ri);
1285                        }
[410]1286                        else
[408]1287                                (*ri).mRay->Unref();
[410]1288                }
1289        }
[465]1290
[410]1291        // update stats
[420]1292        mStat.rayRefs -= (int)leaf->mRays.size();
[414]1293        mStat.rayRefs += raysBack + raysFront;
[408]1294
[420]1295        DEL_PTR(leaf);
1296
[410]1297        return node;
[406]1298}
1299
[410]1300int VspKdTree::ReleaseMemory(const int time)
1301{
[462]1302        stack<VspKdNode *> tstack;
[406]1303
[410]1304        // find a node in the tree which subtree will be collapsed
[422]1305        int maxAccessTime = time - mAccessTimeThreshold;
[420]1306        int released = 0;
[406]1307
[419]1308        tstack.push(mRoot);
[406]1309
[423]1310        while (!tstack.empty())
[410]1311        {
[462]1312                VspKdNode *node = tstack.top();
[410]1313                tstack.pop();
[465]1314
1315                if (!node->IsLeaf())
[410]1316                {
[462]1317                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
[410]1318                        // cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
[465]1319                        if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
[410]1320                        {
1321                                released = CollapseSubtree(node, time);
1322                                break;
1323                        }
[465]1324                        if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
[410]1325                        {
[419]1326                                tstack.push(in->GetFront());
1327                                tstack.push(in->GetBack());
[465]1328                        }
1329                        else
[410]1330                        {
[419]1331                                tstack.push(in->GetBack());
1332                                tstack.push(in->GetFront());
[410]1333                        }
1334                }
1335        }
[406]1336
[465]1337        while (tstack.empty())
[410]1338        {
1339                // could find node to collaps...
1340                // cout<<"Could not find a node to release "<<endl;
1341                break;
1342        }
[465]1343
[410]1344        return released;
[408]1345}
[406]1346
1347
[462]1348VspKdNode *VspKdTree::SubdivideLeaf(VspKdLeaf *leaf,
[410]1349                                                                                const float sizeThreshold)
[408]1350{
[462]1351        VspKdNode *node = leaf;
[406]1352
[410]1353        AxisAlignedBox3 leafBBox = GetBBox(leaf);
[406]1354
[408]1355        static int pass = 0;
[410]1356        ++ pass;
[465]1357
[410]1358        // check if we should perform a dynamic subdivision of the leaf
[420]1359        if (// leaf->mRays.size() > (unsigned)termMinCost &&
[465]1360                (leaf->GetPvsSize() >= mTermMinPvs) &&
[410]1361                (SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
1362        {
[423]1363        // memory check and release
[465]1364                if (GetMemUsage() > mMaxTotalMemory)
[410]1365                        ReleaseMemory(pass);
[465]1366
[410]1367                AxisAlignedBox3 backBBox, frontBBox;
[406]1368
[410]1369                // subdivide the node
1370                node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
1371        }
1372        return node;
[406]1373}
1374
1375
1376
[423]1377void VspKdTree::UpdateRays(VssRayContainer &remove,
1378                                                   VssRayContainer &add)
[406]1379{
[462]1380        VspKdLeaf::NewMail();
[406]1381
[410]1382        // schedule rays for removal
1383        for(VssRayContainer::const_iterator ri = remove.begin();
[465]1384                ri != remove.end(); ++ ri)
[410]1385        {
[465]1386                (*ri)->ScheduleForRemoval();
[410]1387        }
[406]1388
[410]1389        int inactive = 0;
[406]1390
[465]1391        for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
[410]1392        {
1393                if ((*ri)->ScheduledForRemoval())
1394                        //      RemoveRay(*ri, NULL, false);
1395                        // !!! BUG - with true it does not work correctly - aggreated delete
1396                        RemoveRay(*ri, NULL, true);
1397                else
1398                        ++ inactive;
1399        }
[406]1400
[410]1401        //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
[465]1402        for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
[410]1403        {
1404                AddRay(*ri);
1405        }
[406]1406}
1407
1408
[410]1409void VspKdTree::RemoveRay(VssRay *ray,
[462]1410                                                  vector<VspKdLeaf *> *affectedLeaves,
[410]1411                                                  const bool removeAllScheduledRays)
[406]1412{
[410]1413        stack<RayTraversalData> tstack;
[406]1414
[420]1415        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[465]1416
[410]1417        RayTraversalData data;
[406]1418
[410]1419        // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
[465]1420        while (!tstack.empty())
[410]1421        {
1422                data = tstack.top();
1423                tstack.pop();
[406]1424
[465]1425                if (!data.mNode->IsLeaf())
[410]1426                {
1427                        // split the set of rays in two groups intersecting the
1428                        // two subtrees
1429                        TraverseInternalNode(data, tstack);
[465]1430        }
1431                else
[410]1432                {
1433                        // remove the ray from the leaf
1434                        // find the ray in the leaf and swap it with the last ray...
[462]1435                        VspKdLeaf *leaf = (VspKdLeaf *)data.mNode;
[465]1436
1437                        if (!leaf->Mailed())
[410]1438                        {
1439                                leaf->Mail();
[465]1440
[410]1441                                if (affectedLeaves)
1442                                        affectedLeaves->push_back(leaf);
[465]1443
1444                                if (removeAllScheduledRays)
[410]1445                                {
[420]1446                                        int tail = (int)leaf->mRays.size() - 1;
[406]1447
[465]1448                                        for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
[410]1449                                        {
[465]1450                                                if (leaf->mRays[i].mRay->ScheduledForRemoval())
[410]1451                                                {
1452                                                        // find a ray to replace it with
[465]1453                                                        while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
[410]1454                                                        {
[414]1455                                                                ++ mStat.removedRayRefs;
[420]1456                                                                leaf->mRays[tail].mRay->Unref();
1457                                                                leaf->mRays.pop_back();
[465]1458
[410]1459                                                                -- tail;
1460                                                        }
[465]1461
[410]1462                                                        if (tail < i)
1463                                                                break;
[465]1464
[414]1465                                                        ++ mStat.removedRayRefs;
[465]1466
[420]1467                                                        leaf->mRays[i].mRay->Unref();
1468                                                        leaf->mRays[i] = leaf->mRays[tail];
1469                                                        leaf->mRays.pop_back();
[410]1470                                                        -- tail;
1471                                                }
1472                                        }
1473                                }
1474                        }
[465]1475
[410]1476                        if (!removeAllScheduledRays)
[465]1477                                for (int i=0; i < (int)leaf->mRays.size(); i++)
[410]1478                                {
[465]1479                                        if (leaf->mRays[i].mRay == ray)
[410]1480                                        {
[414]1481                                                ++ mStat.removedRayRefs;
[410]1482                                                ray->Unref();
[420]1483                                                leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
1484                                                leaf->mRays.pop_back();
[410]1485                                            // check this ray again
1486                                                break;
1487                                        }
1488                                }
1489                }
[408]1490        }
1491
[465]1492        if (ray->RefCount() != 0)
[410]1493        {
1494                cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1495                exit(1);
1496        }
1497
[406]1498}
1499
[410]1500void VspKdTree::AddRay(VssRay *ray)
[406]1501{
[410]1502        stack<RayTraversalData> tstack;
[465]1503
[420]1504        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[465]1505
[410]1506        RayTraversalData data;
[465]1507
1508        while (!tstack.empty())
[410]1509        {
1510                data = tstack.top();
1511                tstack.pop();
[406]1512
[465]1513                if (!data.mNode->IsLeaf())
[410]1514                {
1515                        TraverseInternalNode(data, tstack);
[465]1516                }
1517                else
[410]1518                {
1519                        // remove the ray from the leaf
[420]1520                        // find the ray in the leaf and swap it with the last ray
[462]1521                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(data.mNode);
[423]1522
[419]1523                        leaf->AddRay(data.mRayData);
[413]1524                        ++ mStat.addedRayRefs;
[410]1525                }
1526        }
[406]1527}
1528
[465]1529void VspKdTree::TraverseInternalNode(RayTraversalData &data,
[410]1530                                                                         stack<RayTraversalData> &tstack)
[406]1531{
[462]1532        VspKdInterior *in = dynamic_cast<VspKdInterior *>(data.mNode);
[465]1533
1534        if (in->mAxis <= VspKdNode::SPLIT_Z)
[410]1535        {
1536            // determine the side of this ray with respect to the plane
[485]1537                float t;
1538                const int side = in->ComputeRayIntersection(data.mRayData, t);
[465]1539
1540                if (side == 0)
[410]1541                {
[465]1542                        if (data.mRayData.mRay->HasPosDir(in->mAxis))
[433]1543                        {
[465]1544                                tstack.push(RayTraversalData(in->GetBack(),
[463]1545                                                        RayInfo(data.mRayData.mRay, data.mRayData.mMinT, t)));
[465]1546
[433]1547                                tstack.push(RayTraversalData(in->GetFront(),
[463]1548                                                        RayInfo(data.mRayData.mRay, t, data.mRayData.mMaxT)));
[465]1549
1550                        }
1551                        else
[433]1552                        {
1553                                tstack.push(RayTraversalData(in->GetBack(),
1554                                                        RayInfo(data.mRayData.mRay,
[463]1555                                                        t,
[433]1556                                                        data.mRayData.mMaxT)));
1557                                tstack.push(RayTraversalData(in->GetFront(),
1558                                                        RayInfo(data.mRayData.mRay,
1559                                                        data.mRayData.mMinT,
[463]1560                                                        t)));
[433]1561                        }
[410]1562                }
[433]1563                else
1564                        if (side == 1)
1565                                tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1566                        else
1567                                tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[465]1568        }
1569        else
[410]1570        {
1571                // directional split
[419]1572                if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
1573                        tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
[408]1574                else
[419]1575                        tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[410]1576        }
[406]1577}
1578
[408]1579
[462]1580int VspKdTree::CollapseSubtree(VspKdNode *sroot, const int time)
[406]1581{
[410]1582        // first count all rays in the subtree
1583        // use mail 1 for this purpose
[462]1584        stack<VspKdNode *> tstack;
[465]1585
[410]1586        int rayCount = 0;
1587        int totalRayCount = 0;
1588        int collapsedNodes = 0;
[408]1589
1590#if DEBUG_COLLAPSE
[420]1591        cout << "Collapsing subtree" << endl;
[428]1592        cout << "accessTime=" << sroot->GetAccessTime() << endl;
[420]1593        cout << "depth=" << (int)sroot->depth << endl;
[408]1594#endif
1595
[410]1596        // tstat.collapsedSubtrees++;
1597        // tstat.collapseDepths += (int)sroot->depth;
1598        // tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1599        tstack.push(sroot);
1600        VssRay::NewMail();
[465]1601
1602        while (!tstack.empty())
[410]1603        {
1604                ++ collapsedNodes;
[465]1605
[462]1606                VspKdNode *node = tstack.top();
[410]1607                tstack.pop();
[465]1608                if (node->IsLeaf())
[410]1609                {
[462]1610                        VspKdLeaf *leaf = (VspKdLeaf *) node;
[465]1611
[420]1612                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1613                                        ri != leaf->mRays.end(); ++ ri)
[410]1614                        {
1615                                ++ totalRayCount;
[465]1616
1617                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed())
[410]1618                                {
[408]1619                                        (*ri).mRay->Mail();
[410]1620                                        ++ rayCount;
[408]1621                                }
[410]1622                        }
[465]1623                }
1624                else
[410]1625                {
[472]1626                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetFront());
1627                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetBack());
[410]1628                }
1629        }
[465]1630
[410]1631        VssRay::NewMail();
[406]1632
[410]1633        // create a new node that will hold the rays
[462]1634        VspKdLeaf *newLeaf = new VspKdLeaf(sroot->mParent, rayCount);
[465]1635
[500]1636        VspKdInterior *parent =
1637                dynamic_cast<VspKdInterior *>(newLeaf->mParent);
1638
1639        if (parent)
[472]1640        {
[500]1641                parent->ReplaceChildLink(sroot, newLeaf);
[472]1642        }
[410]1643        tstack.push(sroot);
[406]1644
[465]1645        while (!tstack.empty())
[410]1646        {
[462]1647                VspKdNode *node = tstack.top();
[410]1648                tstack.pop();
[408]1649
[465]1650                if (node->IsLeaf())
[410]1651                {
[462]1652                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[408]1653
[420]1654                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]1655                                        ri != leaf->mRays.end(); ++ ri)
[410]1656                        {
[465]1657                                // unref this ray from the old node
1658                                if ((*ri).mRay->IsActive())
[410]1659                                {
[408]1660                                        (*ri).mRay->Unref();
[465]1661                                        if (!(*ri).mRay->Mailed())
[410]1662                                        {
[408]1663                                                (*ri).mRay->Mail();
1664                                                newLeaf->AddRay(*ri);
1665                                        }
[465]1666                                }
[410]1667                                else
[408]1668                                        (*ri).mRay->Unref();
[410]1669                        }
[465]1670                }
1671                else
[410]1672                {
[465]1673                        VspKdInterior *interior =
[462]1674                                dynamic_cast<VspKdInterior *>(node);
[419]1675                        tstack.push(interior->GetBack());
1676                        tstack.push(interior->GetFront());
[410]1677                }
1678        }
[465]1679
[410]1680        // delete the node and all its children
[420]1681        DEL_PTR(sroot);
[465]1682
[462]1683        // for(VspKdNode::SRayContainer::iterator ri = newleaf->mRays.begin();
[420]1684    //      ri != newleaf->mRays.end(); ++ ri)
[410]1685        //     (*ri).ray->UnMail(2);
[408]1686
1687#if DEBUG_COLLAPSE
[423]1688        cout << "Total memory before=" << GetMemUsage() << endl;
[408]1689#endif
1690
[414]1691        mStat.nodes -= collapsedNodes - 1;
1692        mStat.rayRefs -= totalRayCount - rayCount;
[465]1693
[408]1694#if DEBUG_COLLAPSE
[410]1695        cout << "collapsed nodes" << collapsedNodes << endl;
1696        cout << "collapsed rays" << totalRayCount - rayCount << endl;
1697        cout << "Total memory after=" << GetMemUsage() << endl;
1698        cout << "================================" << endl;
[408]1699#endif
1700
1701        //  tstat.collapsedNodes += collapsedNodes;
1702        //  tstat.collapsedRays += totalRayCount - rayCount;
[410]1703    return totalRayCount - rayCount;
[406]1704}
1705
[408]1706
[462]1707int VspKdTree::GetPvsSize(VspKdNode *node, const AxisAlignedBox3 &box) const
[406]1708{
[462]1709        stack<VspKdNode *> tstack;
[419]1710        tstack.push(mRoot);
[408]1711
1712        Intersectable::NewMail();
1713        int pvsSize = 0;
[465]1714
1715        while (!tstack.empty())
[410]1716        {
[462]1717                VspKdNode *node = tstack.top();
[410]1718                tstack.pop();
[465]1719
1720          if (node->IsLeaf())
[410]1721                {
[462]1722                        VspKdLeaf *leaf = (VspKdLeaf *)node;
[420]1723                        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1724                                 ri != leaf->mRays.end(); ++ ri)
[410]1725                        {
[465]1726                                if ((*ri).mRay->IsActive())
[410]1727                                {
[408]1728                                        Intersectable *object;
1729#if BIDIRECTIONAL_RAY
1730                                        object = (*ri).mRay->mOriginObject;
[465]1731
1732                                        if (object && !object->Mailed())
[410]1733                                        {
1734                                                ++ pvsSize;
[408]1735                                                object->Mail();
1736                                        }
1737#endif
1738                                        object = (*ri).mRay->mTerminationObject;
[465]1739                                        if (object && !object->Mailed())
[410]1740                                        {
1741                                                ++ pvsSize;
[408]1742                                                object->Mail();
1743                                        }
1744                                }
[465]1745                        }
[410]1746                }
[465]1747                else
[410]1748                {
[462]1749                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
[465]1750
1751                        if (in->mAxis < 3)
[410]1752                        {
[419]1753                                if (box.Max(in->mAxis) >= in->mPosition)
1754                                        tstack.push(in->GetFront());
[465]1755
[419]1756                                if (box.Min(in->mAxis) <= in->mPosition)
1757                                        tstack.push(in->GetBack());
[465]1758                        }
1759                        else
[410]1760                        {
[408]1761                                // both nodes for directional splits
[419]1762                                tstack.push(in->GetFront());
1763                                tstack.push(in->GetBack());
[408]1764                        }
1765                }
1766        }
[410]1767
[408]1768        return pvsSize;
[406]1769}
1770
[410]1771void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1772                                                                                         float &maxRayContribution,
1773                                                                                         float &avgRayContribution)
[406]1774{
[462]1775        stack<VspKdNode *> tstack;
[419]1776        tstack.push(mRoot);
[465]1777
[408]1778        minRayContribution = 1.0f;
1779        maxRayContribution = 0.0f;
[410]1780
[408]1781        float sumRayContribution = 0.0f;
1782        int leaves = 0;
[406]1783
[465]1784        while (!tstack.empty())
[410]1785        {
[462]1786                VspKdNode *node = tstack.top();
[410]1787                tstack.pop();
[465]1788                if (node->IsLeaf())
[410]1789                {
[408]1790                        leaves++;
[462]1791                        VspKdLeaf *leaf = (VspKdLeaf *)node;
[408]1792                        float c = leaf->GetAvgRayContribution();
[410]1793
[408]1794                        if (c > maxRayContribution)
1795                                maxRayContribution = c;
1796                        if (c < minRayContribution)
1797                                minRayContribution = c;
1798                        sumRayContribution += c;
[465]1799
1800                }
[482]1801                else
1802                {
[462]1803                        VspKdInterior *in = (VspKdInterior *)node;
[482]1804                       
[419]1805                        tstack.push(in->GetFront());
1806                        tstack.push(in->GetBack());
[408]1807                }
1808        }
[465]1809
[423]1810        Debug << "sum=" << sumRayContribution << endl;
1811        Debug << "leaves=" << leaves << endl;
1812        avgRayContribution = sumRayContribution / (float)leaves;
[406]1813}
1814
1815
[410]1816int VspKdTree::GenerateRays(const float ratioPerLeaf,
1817                                                        SimpleRayContainer &rays)
[406]1818{
[462]1819        stack<VspKdNode *> tstack;
[419]1820        tstack.push(mRoot);
[465]1821
1822        while (!tstack.empty())
[410]1823        {
[462]1824                VspKdNode *node = tstack.top();
[410]1825                tstack.pop();
[465]1826
1827                if (node->IsLeaf())
[410]1828                {
[462]1829                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[423]1830
[408]1831                        float c = leaf->GetAvgRayContribution();
[411]1832                        int num = (int)(c*ratioPerLeaf + 0.5);
[408]1833                        //                      cout<<num<<" ";
1834
[465]1835                        for (int i=0; i < num; i++)
[410]1836                        {
[408]1837                                Vector3 origin = GetBBox(leaf).GetRandomPoint();
[419]1838                                /*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
[408]1839                                Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1840                                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[419]1841                                rays.push_back(SimpleRay(origin, direction));*/
[408]1842                        }
[465]1843
1844                }
1845                else
[410]1846                {
[465]1847                        VspKdInterior *in =
[462]1848                                dynamic_cast<VspKdInterior *>(node);
[482]1849               
[419]1850                        tstack.push(in->GetFront());
1851                        tstack.push(in->GetBack());
[408]1852                }
[406]1853        }
1854
[411]1855        return (int)rays.size();
[406]1856}
1857
1858
[410]1859float VspKdTree::GetAvgPvsSize()
[406]1860{
[462]1861        stack<VspKdNode *> tstack;
[419]1862        tstack.push(mRoot);
[406]1863
[408]1864        int sumPvs = 0;
1865        int leaves = 0;
[465]1866
1867        while (!tstack.empty())
[410]1868        {
[462]1869                VspKdNode *node = tstack.top();
[410]1870                tstack.pop();
[465]1871
1872                if (node->IsLeaf())
[410]1873                {
[462]1874                        VspKdLeaf *leaf = (VspKdLeaf *)node;
[408]1875                        // update pvs size
1876                        leaf->UpdatePvsSize();
1877                        sumPvs += leaf->GetPvsSize();
1878                        leaves++;
[465]1879                }
1880                else
[410]1881                {
[462]1882                        VspKdInterior *in = (VspKdInterior *)node;
[482]1883                       
[419]1884                        tstack.push(in->GetFront());
1885                        tstack.push(in->GetBack());
[406]1886                }
1887        }
[408]1888
[410]1889        return sumPvs / (float)leaves;
[418]1890}
1891
[462]1892VspKdNode *VspKdTree::GetRoot() const
[418]1893{
1894        return mRoot;
1895}
1896
[462]1897AxisAlignedBox3 VspKdTree::GetBBox(VspKdNode *node) const
[418]1898{
[419]1899        if (node->mParent == NULL)
1900                return mBox;
[418]1901
1902        if (!node->IsLeaf())
[500]1903        {
1904                if (node->Type() == VspKdNode::EIntermediate)
1905            return (dynamic_cast<VspKdInterior *>(node))->mBox;
1906                else
1907                        return (dynamic_cast<VspKdIntermediate *>(node))->mBox;
1908        }
[418]1909
[500]1910        VspKdInterior *parent = dynamic_cast<VspKdInterior *>(node->mParent);
[465]1911
[500]1912        if (parent->mAxis >= 3)
1913                return parent->mBox;
1914
1915        AxisAlignedBox3 box(parent->mBox);
1916        if (parent->GetFront() == node)
1917                box.SetMin(parent->mAxis, parent->mPosition);
[418]1918        else
[500]1919                box.SetMax(parent->mAxis, parent->mPosition);
[419]1920
1921        return box;
[418]1922}
1923
[465]1924int     VspKdTree::GetRootPvsSize() const
[418]1925{
1926        return GetPvsSize(mRoot, mBox);
1927}
[419]1928
[465]1929const VspKdStatistics &VspKdTree::GetStatistics() const
[419]1930{
1931        return mStat;
[420]1932}
1933
1934void VspKdTree::AddRays(VssRayContainer &add)
1935{
1936        VssRayContainer remove;
1937        UpdateRays(remove, add);
1938}
[465]1939
[420]1940// return memory usage in MB
[465]1941float VspKdTree::GetMemUsage() const
[420]1942{
[465]1943        return
1944                (sizeof(VspKdTree) +
[462]1945                 mStat.Leaves() * sizeof(VspKdLeaf) +
1946                 mStat.Interior() * sizeof(VspKdInterior) +
[420]1947                 mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
1948}
[465]1949
1950float VspKdTree::GetRayMemUsage() const
[420]1951{
1952        return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
[425]1953}
1954
[482]1955
[465]1956int VspKdTree::ComputePvsSize(VspKdNode *node,
[426]1957                                                          const RayInfoContainer &globalRays) const
[425]1958{
1959        int pvsSize = 0;
1960
1961        const AxisAlignedBox3 box = GetBBox(node);
1962
1963        RayInfoContainer::const_iterator it, it_end = globalRays.end();
1964
[495]1965        Intersectable::NewMail();
[482]1966
[426]1967        // warning: implicit conversion from VssRay to Ray
[425]1968        for (it = globalRays.begin(); it != globalRays.end(); ++ it)
[495]1969        {
1970                VssRay *ray = (*it).mRay;
[465]1971
[495]1972                // ray intersects node bounding box
1973                if (box.GetMinMaxT(*ray, NULL, NULL))
1974                {
1975                        if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
1976                        {
1977                                ray->mTerminationObject->Mail();
1978                                ++ pvsSize;
1979                        }
1980                }
1981        }
[425]1982        return pvsSize;
[428]1983}
1984
[482]1985
[462]1986void VspKdTree::CollectLeaves(vector<VspKdLeaf *> &leaves) const
[428]1987{
[462]1988        stack<VspKdNode *> nodeStack;
[428]1989        nodeStack.push(mRoot);
1990
[465]1991        while (!nodeStack.empty())
[428]1992        {
[462]1993                VspKdNode *node = nodeStack.top();
[465]1994
[428]1995                nodeStack.pop();
[465]1996
1997                if (node->IsLeaf())
[428]1998                {
[462]1999                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[428]2000                        leaves.push_back(leaf);
[465]2001                }
2002                else
[428]2003                {
[462]2004                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
[428]2005
[448]2006                        nodeStack.push(interior->GetBack());
2007                        nodeStack.push(interior->GetFront());
[428]2008                }
2009        }
[452]2010}
2011
[503]2012
[462]2013int VspKdTree::FindNeighbors(VspKdLeaf *n,
2014                                                         vector<VspKdLeaf *> &neighbors,
[452]2015                                                         bool onlyUnmailed)
2016{
[462]2017        stack<VspKdNode *> nodeStack;
[465]2018
[452]2019        nodeStack.push(mRoot);
2020
2021        AxisAlignedBox3 box = GetBBox(n);
2022
[465]2023        while (!nodeStack.empty())
[452]2024        {
[462]2025                VspKdNode *node = nodeStack.top();
[452]2026                nodeStack.pop();
2027
[465]2028                if (node->IsLeaf())
[452]2029                {
[462]2030                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
[452]2031
[465]2032                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
[452]2033                                neighbors.push_back(leaf);
[465]2034                }
2035                else
[452]2036                {
[462]2037                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
[452]2038
2039                        if (interior->mPosition > box.Max(interior->mAxis))
2040                                nodeStack.push(interior->mBack);
2041                        else
2042                        {
2043                                if (interior->mPosition < box.Min(interior->mAxis))
2044                                        nodeStack.push(interior->mFront);
[465]2045                                else
[452]2046                                {
2047                                        // random decision
2048                                        nodeStack.push(interior->mBack);
2049                                        nodeStack.push(interior->mFront);
2050                                }
2051                        }
2052                }
2053        }
2054
2055        return (int)neighbors.size();
2056}
2057
[462]2058
2059void VspKdTree::SetViewCellsManager(ViewCellsManager *vcm)
[452]2060{
[462]2061        mViewCellsManager = vcm;
2062}
[452]2063
[462]2064
[468]2065int VspKdTree::CastLineSegment(const Vector3 &origin,
2066                                                           const Vector3 &termination,
2067                                                           vector<ViewCell *> &viewcells)
[462]2068{
2069        int hits = 0;
[465]2070
[468]2071        float mint = 0.0f, maxt = 1.0f;
2072        const Vector3 dir = termination - origin;
[465]2073
[483]2074        stack<LineTraversalData> tStack;
[462]2075
2076        Intersectable::NewMail();
2077
[468]2078        Vector3 entp = origin;
2079        Vector3 extp = termination;
[465]2080
[462]2081        VspKdNode *node = mRoot;
2082        VspKdNode *farChild;
2083
2084        float position;
2085        int axis;
[465]2086
2087        while (1)
[462]2088        {
[465]2089                if (!node->IsLeaf())
[462]2090                {
2091                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
2092                        position = in->mPosition;
2093                        axis = in->mAxis;
2094
[465]2095                        if (entp[axis] <= position)
[462]2096                        {
[465]2097                                if (extp[axis] <= position)
[462]2098                                {
2099                                        node = in->mBack;
2100                                        // cases N1,N2,N3,P5,Z2,Z3
2101                                        continue;
[465]2102                                } else
[462]2103                                {
2104                                        // case N4
2105                                        node = in->mBack;
2106                                        farChild = in->mFront;
2107                                }
2108                        }
[465]2109                        else
[462]2110                        {
[465]2111                                if (position <= extp[axis])
[462]2112                                {
2113                                        node = in->mFront;
2114                                        // cases P1,P2,P3,N5,Z1
2115                                        continue;
[465]2116                                }
2117                                else
[462]2118                                {
2119                                        node = in->mFront;
2120                                        farChild = in->mBack;
2121                                        // case P4
2122                                }
2123                        }
2124
2125                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2126                        // case N4 or P4
[468]2127                        float tdist = (position - origin[axis]) / dir[axis];
[483]2128                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
[468]2129                        extp = origin + dir * tdist;
[462]2130                        maxt = tdist;
[465]2131                }
2132                else
[462]2133                {
2134                        // compute intersection with all objects in this leaf
[468]2135                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2136                        ViewCell *vc = leaf->GetViewCell();
[465]2137
[468]2138                        if (!vc->Mailed())
[462]2139                        {
[468]2140                                vc->Mail();
2141                                viewcells.push_back(vc);
[469]2142                                ++ hits;
[462]2143                        }
[465]2144
[485]2145                        if (0) //matt TODO: REMOVE LATER
[483]2146                                leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2147
[462]2148                        // get the next node from the stack
2149                        if (tStack.empty())
2150                                break;
[465]2151
[462]2152                        entp = extp;
2153                        mint = maxt;
[468]2154                       
[483]2155                        LineTraversalData &s  = tStack.top();
[462]2156                        node = s.mNode;
2157                        extp = s.mExitPoint;
2158                        maxt = s.mMaxT;
2159                        tStack.pop();
2160                }
2161        }
[465]2162
[462]2163        return hits;
2164}
2165
2166
2167void VspKdTree::GenerateViewCell(VspKdLeaf *leaf)
2168{
2169        RayInfoContainer::const_iterator it, it_end = leaf->GetRays().end();
2170
2171        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(mViewCellsManager->GenerateViewCell());
2172        leaf->SetViewCell(vc);
2173
[479]2174        vc->SetVolume(GetBBox(leaf).GetVolume());
2175        vc->SetArea(GetBBox(leaf).SurfaceArea());
2176
[580]2177        vc->mLeaf = leaf;
[462]2178
2179        for (it = leaf->GetRays().begin(); it != it_end; ++ it)
2180        {
2181                VssRay *ray = (*it).mRay;
2182
[463]2183                if (ray->mTerminationObject)
[556]2184                  vc->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf);
[465]2185
[463]2186                if (ray->mOriginObject)
[556]2187                  vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf);
[462]2188        }
2189}
2190
2191
[589]2192bool VspKdTree::MergeViewCells(VspKdLeaf *l1, VspKdLeaf *l2)
2193{
2194        //-- change pointer to view cells of all leaves associated
2195        //-- with the previous view cells
2196        VspKdViewCell *fVc = l1->GetViewCell();
2197        VspKdViewCell *bVc = l2->GetViewCell();
2198
2199        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(
2200                mViewCellsManager->MergeViewCells(fVc, bVc));
2201
2202        // if merge was unsuccessful
2203        if (!vc) return false;
2204// TODO
2205#if VC_HISTORY
2206        // set new size of view cell
[479]2207        vc->SetVolume(fVc->GetVolume() + bVc->GetVolume());
[482]2208        // new area
[589]2209        vc->SetArea(fVc->GetArea() + bVc->GetArea());
2210
2211        vector<VspKdLeaf *> fLeaves = fVc->mLeaves;
2212        vector<VspKdLeaf *> bLeaves = bVc->mLeaves;
2213
2214        vector<VspKdLeaf *>::const_iterator it;
2215
2216        //-- change view cells of all the other leaves the view cell belongs to
2217        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
2218        {
2219                (*it)->SetViewCell(vc);
2220                vc->mLeaves.push_back(*it);
2221        }
2222
2223        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
2224        {
2225                (*it)->SetViewCell(vc);
2226                vc->mLeaves.push_back(*it);
2227        }
2228#endif
[465]2229        vc->Mail();
[589]2230
2231        // clean up old view cells
2232        DEL_PTR(fVc);
2233        DEL_PTR(bVc);
2234
[463]2235        return true;
[462]2236}
2237
[483]2238void VspKdTree::CollectMergeCandidates()
[462]2239{
[580]2240        VspKdMergeCandidate::sOverallCost = 0;
[462]2241        vector<VspKdLeaf *> leaves;
[465]2242
[462]2243        // collect the leaves, e.g., the "voxels" that will build the view cells
[452]2244        CollectLeaves(leaves);
2245
[462]2246        VspKdLeaf::NewMail();
[452]2247
[462]2248        vector<VspKdLeaf *>::const_iterator it, it_end = leaves.end();
[452]2249
[462]2250        // generate view cells
[452]2251        for (it = leaves.begin(); it != it_end; ++ it)
[462]2252                GenerateViewCell(*it);
2253
2254        // find merge candidates and push them into queue
2255        for (it = leaves.begin(); it != it_end; ++ it)
[452]2256        {
[462]2257                VspKdLeaf *leaf = *it;
[485]2258                ViewCell *vc = leaf->GetViewCell();
[462]2259                // no leaf is part of two merge candidates
[485]2260                leaf->Mail();
[580]2261                VspKdMergeCandidate::sOverallCost +=
[485]2262                        vc->GetVolume() * vc->GetPvs().GetSize();
2263                vector<VspKdLeaf *> neighbors;
2264                FindNeighbors(leaf, neighbors, true);
[462]2265
[485]2266                vector<VspKdLeaf *>::const_iterator nit,
2267                        nit_end = neighbors.end();
[452]2268
[485]2269                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
2270                {
[580]2271                        mMergeQueue.push(VspKdMergeCandidate(leaf, *nit));
[462]2272                }
2273        }
[483]2274}
2275
2276
[485]2277int VspKdTree::MergeViewCells(const VssRayContainer &rays)
[483]2278{
[485]2279        CollectMergeCandidates();
2280
[462]2281        int merged = 0;
2282
[485]2283        int vcSize = mStat.Leaves();
[462]2284        // use priority queue to merge leaves
[483]2285        while (!mMergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
2286                   (mMergeQueue.top().GetMergeCost() <
[580]2287                    mMergeMaxCostRatio * VspKdMergeCandidate::sOverallCost))
[462]2288        {
[485]2289                //Debug << "mergecost: " << mergeQueue.top().GetMergeCost() /
[580]2290                //VspKdMergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
2291                VspKdMergeCandidate mc = mMergeQueue.top();
[483]2292                mMergeQueue.pop();
[465]2293
2294                // both view cells equal!
[462]2295                if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
2296                        continue;
2297
2298                if (mc.Valid())
2299                {
[465]2300                        ViewCell::NewMail();
2301                        MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
[462]2302
2303                        ++ merged;
2304                        -- vcSize;
[472]2305                        // increase absolute merge cost
[580]2306                        VspKdMergeCandidate::sOverallCost += mc.GetMergeCost();
[452]2307                }
[462]2308                // merge candidate not valid, because one of the leaves was already
2309                // merged with another one
[465]2310                else
[462]2311                {
2312                        // validate and reinsert into queue
2313                        mc.SetValid();
[483]2314                        mMergeQueue.push(mc);
[465]2315                        //Debug << "validate " << mc.GetMergeCost() << endl;
[462]2316                }
[452]2317        }
[462]2318
[485]2319        Debug << "merged " << merged << " of " << mStat.Leaves() << " leaves" << endl;
[465]2320
[469]2321        //TODO: should return sample contributions
[462]2322        return merged;
[452]2323}
2324
2325
[517]2326void VspKdTree::RepairViewCellsLeafLists()
[589]2327{
[465]2328        // list not valid anymore => clear
[589]2329        stack<VspKdNode *> nodeStack;
2330        nodeStack.push(mRoot);
2331
2332        ViewCell::NewMail();
2333
2334        while (!nodeStack.empty())
2335        {
2336                VspKdNode *node = nodeStack.top();
2337                nodeStack.pop();
2338
2339                if (node->IsLeaf())
2340                {
2341                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2342
2343                        VspKdViewCell *viewCell = leaf->GetViewCell();
2344#if VC_HISTORY
2345                        if (!viewCell->Mailed())
2346                        {
2347                                viewCell->mLeaves.clear();
2348                                viewCell->Mail();
2349                        }
2350
2351                        viewCell->mLeaves.push_back(leaf);
2352#endif
2353                }
2354                else
2355                {
2356                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2357
2358                        nodeStack.push(interior->GetFront());
2359                        nodeStack.push(interior->GetBack());
2360                }
[465]2361        }
2362}
2363
[479]2364
[501]2365VspKdNode *VspKdTree::CollapseTree(VspKdNode *node, int &collapsed)
[465]2366{
[501]2367        if (node->IsLeaf())
[465]2368                return node;
2369
[501]2370        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2371
2372        VspKdNode *front = CollapseTree(interior->GetFront(), collapsed);
2373        VspKdNode *back = CollapseTree(interior->GetBack(), collapsed);
2374
[465]2375        if (front->IsLeaf() && back->IsLeaf())
2376        {
2377                VspKdLeaf *frontLeaf = dynamic_cast<VspKdLeaf *>(front);
2378                VspKdLeaf *backLeaf = dynamic_cast<VspKdLeaf *>(back);
2379
2380                //-- collapse tree
2381                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
2382                {
2383                        VspKdViewCell *vc = frontLeaf->GetViewCell();
2384
2385                        VspKdLeaf *leaf = new VspKdLeaf(interior->GetParent(), 0);
2386                        leaf->SetViewCell(vc);
[501]2387                        //leaf->SetTreeValid(frontLeaf->TreeValid());
[465]2388
2389                        // replace a link from node's parent
[501]2390                        if (leaf->mParent)
2391                        {
2392                                // replace a link from node's parent
2393                                VspKdInterior *parent =
2394                                        dynamic_cast<VspKdInterior *>(leaf->mParent);
[500]2395                                parent->ReplaceChildLink(node, leaf);
[501]2396                        }
2397                        ++ collapsed;
[465]2398                        delete interior;
2399
2400                        return leaf;
2401                }
2402        }
2403
[495]2404        return node;
2405}
2406
2407
[501]2408int VspKdTree::CollapseTree()
[495]2409{
[501]2410        int collapsed = 0;
2411        CollapseTree(mRoot, collapsed);
[485]2412        // revalidate leaves
[517]2413        RepairViewCellsLeafLists();
[501]2414
2415        return collapsed;
[465]2416}
2417
[501]2418
[486]2419int VspKdTree::RefineViewCells(const VssRayContainer &rays)
[473]2420{
[486]2421        int shuffled = 0;
2422
2423        Debug << "refining " << (int)mMergeQueue.size() << " candidates " << endl;
2424        BspLeaf::NewMail();
2425
2426        // Use priority queue of remaining leaf pairs
2427        // These candidates either share the same view cells or
2428        // are border leaves which share a boundary.
2429        // We test if they can be shuffled, i.e.,
2430        // either one leaf is made part of one view cell or the other
2431        // leaf is made part of the other view cell. It is tested if the
2432        // remaining view cells are "better" than the old ones.
2433        while (!mMergeQueue.empty())
2434        {
[580]2435                VspKdMergeCandidate mc = mMergeQueue.top();
[486]2436                mMergeQueue.pop();
2437
2438                // both view cells equal or already shuffled
2439                if ((mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell()) ||
2440                        (mc.GetLeaf1()->Mailed()) || (mc.GetLeaf2()->Mailed()))
2441                        continue;
2442               
2443                // candidate for shuffling
2444                const bool wasShuffled = false;
2445                        //ShuffleLeaves(mc.GetLeaf1(), mc.GetLeaf2());
[495]2446                        // matt: restore
[486]2447                //-- stats
2448                if (wasShuffled)
2449                        ++ shuffled;
2450        }
2451
2452        return shuffled;
[473]2453}
2454
2455
[486]2456inline int AddedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2457{
2458        return pvs1.AddPvs(pvs2);
2459}
2460
2461
2462inline int SubtractedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2463{
2464        return pvs1.SubtractPvs(pvs2);
2465}
2466
2467
[580]2468float EvalShuffleCost(VspKdLeaf *leaf, VspKdViewCell *vc1, VspKdViewCell *vc2)
[486]2469{
2470        const int pvs1 = SubtractedPvsSize(vc1->GetPvs(), *leaf->mPvs);
2471        const int pvs2 = AddedPvsSize(vc2->GetPvs(), *leaf->mPvs);
2472
[495]2473        const float vol1 = vc1->GetVolume() - leaf->mVolume;
2474        const float vol2 = vc2->GetVolume() + leaf->mVolume;
[486]2475
[495]2476        const float cost1 = pvs1 * vol1;
2477        const float cost2 = pvs2 * vol2;
[486]2478
2479        return cost1 + cost2;
2480}
2481
2482
2483void VspKdTree::ShuffleLeaf(VspKdLeaf *leaf,
2484                                                        VspKdViewCell *vc1,
2485                                                        VspKdViewCell *vc2) const
2486{
2487        // compute new pvs and area
[580]2488        //      TODO
2489#if VC_HISTORY
[486]2490        vc1->GetPvs().SubtractPvs(*leaf->mPvs);
2491        vc2->GetPvs().AddPvs(*leaf->mPvs);
2492       
[495]2493        vc1->SetArea(vc1->GetVolume() - leaf->mVolume);
2494        vc2->SetArea(vc2->GetVolume() + leaf->mVolume);
[486]2495
2496        /// add to second view cell
2497        vc2->mLeaves.push_back(leaf);
2498
2499        // erase leaf from old view cell
2500        vector<VspKdLeaf *>::iterator it = vc1->mLeaves.begin();
2501
2502        for (; *it != leaf; ++ it);
2503        vc1->mLeaves.erase(it);
2504
2505        leaf->SetViewCell(vc2);  // finally change view cell
[580]2506#endif
[486]2507}
2508
2509
2510bool VspKdTree::ShuffleLeaves(VspKdLeaf *leaf1, VspKdLeaf *leaf2) const
2511{
2512        VspKdViewCell *vc1 = leaf1->GetViewCell();
2513        VspKdViewCell *vc2 = leaf2->GetViewCell();
2514
2515        const float cost1 = vc1->GetPvs().GetSize() * vc1->GetArea();
2516        const float cost2 = vc2->GetPvs().GetSize() * vc2->GetArea();
2517
2518        const float oldCost = cost1 + cost2;
2519       
2520        float shuffledCost1 = Limits::Infinity;
2521        float shuffledCost2 = Limits::Infinity;
2522
2523        // the view cell should not be empty after the shuffle
[580]2524//todo
2525        /*      if (vc1->mLeaves.size() > 1)
2526                shuffledCost1 = EvalShuffleCost(leaf1, vc1, vc2);
[486]2527        if (vc2->mLeaves.size() > 1)
[580]2528                shuffledCost2 = EvalShuffleCost(leaf2, vc2, vc1);
[486]2529
[580]2530*/      // shuffling unsuccessful
[486]2531        if ((oldCost <= shuffledCost1) && (oldCost <= shuffledCost2))
2532                return false;
2533       
2534        if (shuffledCost1 < shuffledCost2)
2535        {
2536                ShuffleLeaf(leaf1, vc1, vc2);
2537                leaf1->Mail();
2538        }
2539        else
2540        {
2541                ShuffleLeaf(leaf2, vc2, vc1);
2542                leaf2->Mail();
2543        }
2544
2545        return true;
2546}
2547
2548
2549
2550
[452]2551/*********************************************************************/
[580]2552/*                VspKdMergeCandidate implementation                      */
[452]2553/*********************************************************************/
2554
2555
[580]2556VspKdMergeCandidate::VspKdMergeCandidate(VspKdLeaf *l1, VspKdLeaf *l2):
[462]2557mMergeCost(0),
2558mLeaf1(l1),
2559mLeaf2(l2),
2560mLeaf1Id(l1->GetViewCell()->mMailbox),
2561mLeaf2Id(l2->GetViewCell()->mMailbox)
2562{
2563        EvalMergeCost();
2564}
[452]2565
[580]2566float VspKdMergeCandidate::GetLeaf1Cost() const
[472]2567{
2568        VspKdViewCell *vc = mLeaf1->GetViewCell();
2569        return vc->GetPvs().GetSize() * vc->GetVolume();
2570}
[462]2571
[580]2572float VspKdMergeCandidate::GetLeaf2Cost() const
[472]2573{
2574        VspKdViewCell *vc = mLeaf2->GetViewCell();
2575        return vc->GetPvs().GetSize() * vc->GetVolume();
2576}
2577
[580]2578void VspKdMergeCandidate::EvalMergeCost()
[452]2579{
[465]2580        //-- compute pvs difference
[463]2581        VspKdViewCell *vc1 = mLeaf1->GetViewCell();
[462]2582        VspKdViewCell *vc2 = mLeaf2->GetViewCell();
[452]2583
[463]2584        const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
[462]2585        const int vcPvs = diff1 + vc1->GetPvs().GetSize();
[453]2586
[465]2587        //-- compute ratio of old cost
[462]2588        //-- (added size of left and right view cell times pvs size)
2589        //-- to new rendering cost (size of merged view cell times pvs size)
[472]2590        const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();
2591       
[465]2592        const float newCost =
[469]2593                (float)vcPvs * (vc1->GetVolume() + vc2->GetVolume());
[453]2594
[465]2595        mMergeCost = newCost - oldCost;
[462]2596
[465]2597//      if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
2598//              mMergeCost += 1.0;
[462]2599}
2600
[580]2601void VspKdMergeCandidate::SetLeaf1(VspKdLeaf *l)
[462]2602{
2603        mLeaf1 = l;
2604}
2605
[580]2606void VspKdMergeCandidate::SetLeaf2(VspKdLeaf *l)
[462]2607{
2608        mLeaf2 = l;
2609}
2610
[482]2611
[580]2612VspKdLeaf *VspKdMergeCandidate::GetLeaf1()
[462]2613{
2614        return mLeaf1;
2615}
2616
[482]2617
[580]2618VspKdLeaf *VspKdMergeCandidate::GetLeaf2()
[462]2619{
2620        return mLeaf2;
2621}
2622
[482]2623
[580]2624bool VspKdMergeCandidate::Valid() const
[465]2625{
2626        return
2627                (mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
[462]2628                (mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
2629}
2630
[482]2631
[580]2632float VspKdMergeCandidate::GetMergeCost() const
[462]2633{
2634        return mMergeCost;
2635}
2636
[482]2637
[580]2638void VspKdMergeCandidate::SetValid()
[462]2639{
2640        mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
2641        mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;
[465]2642
[462]2643        EvalMergeCost();
[463]2644}
Note: See TracBrowser for help on using the repository browser.