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

Revision 2015, 59.8 KB checked in by bittner, 17 years ago (diff)

pvs efficiency tuning

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