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

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