source: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp @ 486

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