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

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

pvs efficiency tuning

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