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

Revision 491, 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
834
835float VspKdTree::BestCostRatioRegular(VspKdLeaf *leaf,
836                                                                          const AxisAlignedBox3 &box,
837                                                                          int &axis,
838                                                                          float &position,
839                                                                          int &raysBack,
840                                                                          int &raysFront,
841                                                                          int &pvsBack,
842                                                                          int &pvsFront)
843{
844        int nRaysBack[3], nRaysFront[3];
845        int nPvsBack[3], nPvsFront[3];
846
847        float nPosition[3];
848        float nCostRatio[3];
849        int bestAxis = -1;
850
851        //AxisAlignedBox3 sBox = GetBBox(leaf);
852        const int sAxis = box.Size().DrivingAxis();
853
854        for (axis = 0; axis < 3; ++ axis)
855        {
856                if (!mOnlyDrivingAxis || axis == sAxis)
857                {
858
859                        nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
860
861                        nCostRatio[axis] = EvalCostRatio(leaf,
862                                                                                         box,
863                                                                                         axis,
864                                                                                         nPosition[axis],
865                                                                                         nRaysBack[axis],
866                                                                                         nRaysFront[axis],
867                                                                                         nPvsBack[axis],
868                                                                                         nPvsFront[axis]);
869                        if (bestAxis == -1)
870                        {
871                                bestAxis = axis;
872                        }
873                        else if (nCostRatio[axis] < nCostRatio[bestAxis])
874                        {
875                                /*Debug << "pvs front " << nPvsBack[axis]
876                                          << " pvs back " << nPvsFront[axis]
877                                          << " overall pvs " << leaf->GetPvsSize() << endl;*/
878                                bestAxis = axis;
879                        }
880
881                }
882        }
883
884        //-- assign best axis
885        axis = bestAxis;
886        position = nPosition[bestAxis];
887
888        raysBack = nRaysBack[bestAxis];
889        raysFront = nRaysFront[bestAxis];
890
891        pvsBack = nPvsBack[bestAxis];
892        pvsFront = nPvsFront[bestAxis];
893
894        return nCostRatio[bestAxis];
895}
896
897float VspKdTree::BestCostRatioHeuristic(VspKdLeaf *leaf,
898                                                                                const AxisAlignedBox3 &box,
899                                                                                int &axis,
900                                                                                float &position,
901                                                                                int &raysBack,
902                                                                                int &raysFront,
903                                                                                int &pvsBack,
904                                                                                int &pvsFront)
905{
906        //AxisAlignedBox3 box = GetBBox(leaf);
907        //      AxisAlignedBox3 dirBox = GetDirBBox(node);
908
909        axis = box.Size().DrivingAxis();
910
911        SortSplitCandidates(leaf, axis);
912
913        // go through the lists, count the number of objects left and right
914        // and evaluate the following cost funcion:
915        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
916
917        int rl = 0, rr = (int)leaf->mRays.size();
918        int pl = 0, pr = leaf->GetPvsSize();
919
920        float minBox = box.Min(axis);
921        float maxBox = box.Max(axis);
922        float sizeBox = maxBox - minBox;
923
924        float minBand = minBox + 0.1f*(maxBox - minBox);
925        float maxBand = minBox + 0.9f*(maxBox - minBox);
926
927        float sum = rr*sizeBox;
928        float minSum = 1e20f;
929
930        Intersectable::NewMail();
931
932        // set all object as belonging to the fron pvs
933        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
934                ri != leaf->mRays.end(); ++ ri)
935        {
936                if ((*ri).mRay->IsActive())
937                {
938                        Intersectable *object = (*ri).mRay->mTerminationObject;
939
940                        if (object)
941                                if (!object->Mailed())
942                                {
943                                        object->Mail();
944                                        object->mCounter = 1;
945                                }
946                                else
947                                        ++ object->mCounter;
948                }
949        }
950
951        Intersectable::NewMail();
952
953        for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
954                ci < mSplitCandidates->end(); ++ ci)
955        {
956                VssRay *ray;
957
958                switch ((*ci).type)
959                {
960                        case SortableEntry::ERayMin:
961                                {
962                                        ++ rl;
963                                        ray = (*ci).ray;
964                                        Intersectable *object = ray->mTerminationObject;
965                                        if (object && !object->Mailed())
966                                        {
967                                                object->Mail();
968                                                ++ pl;
969                                        }
970                                        break;
971                                }
972                        case SortableEntry::ERayMax:
973                                {
974                                        -- rr;
975
976                                        ray = (*ci).ray;
977                                        Intersectable *object = ray->mTerminationObject;
978
979                                        if (object)
980                                        {
981                                                if (-- object->mCounter == 0)
982                                                -- pr;
983                                        }
984
985                                        break;
986                                }
987                }
988
989                if ((*ci).value > minBand && (*ci).value < maxBand)
990                {
991                        sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);
992
993                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
994                        // cout<<"cost= "<<sum<<endl;
995
996                        if (sum < minSum)
997                        {
998                                minSum = sum;
999                                position = (*ci).value;
1000
1001                                raysBack = rl;
1002                                raysFront = rr;
1003
1004                                pvsBack = pl;
1005                                pvsFront = pr;
1006
1007                        }
1008                }
1009        }
1010
1011        float oldCost = (float)leaf->GetPvsSize();
1012        float newCost = mCtDivCi + minSum / sizeBox;
1013        float ratio = newCost / oldCost;
1014
1015        //Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1016        //     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
1017
1018        return ratio;
1019}
1020
1021void VspKdTree::SortSplitCandidates(VspKdLeaf *node,
1022                                                                        const int axis)
1023{
1024        mSplitCandidates->clear();
1025
1026        int requestedSize = 2 * (int)(node->mRays.size());
1027        // creates a sorted split candidates array
1028        if (mSplitCandidates->capacity() > 500000 &&
1029                requestedSize < (int)(mSplitCandidates->capacity()/10) )
1030        {
1031        DEL_PTR(mSplitCandidates);
1032                mSplitCandidates = new vector<SortableEntry>;
1033        }
1034
1035        mSplitCandidates->reserve(requestedSize);
1036
1037        // insert all queries
1038        for(RayInfoContainer::const_iterator ri = node->mRays.begin();
1039                ri < node->mRays.end(); ++ ri)
1040        {
1041                bool positive = (*ri).mRay->HasPosDir(axis);
1042                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
1043                                                                                                  (*ri).ExtrapOrigin(axis), (*ri).mRay));
1044
1045                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
1046                                                                                                  (*ri).ExtrapTermination(axis), (*ri).mRay));
1047        }
1048
1049        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
1050}
1051
1052
1053void VspKdTree::EvaluateLeafStats(const TraversalData &data)
1054{
1055        // the node became a leaf -> evaluate stats for leafs
1056        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(data.mNode);
1057
1058        if (leaf->GetPvsSize() > mStat.maxPvsSize)
1059                mStat.maxPvsSize = leaf->GetPvsSize();
1060
1061        if ((int)(leaf->mRays.size()) > mStat.maxRayRefs)
1062                mStat.maxRayRefs = (int)leaf->mRays.size();
1063
1064        mStat.rays += (int)leaf->mRays.size();
1065
1066        if (data.mDepth >= mTermMaxDepth)
1067                ++ mStat.maxDepthNodes;
1068
1069        if (leaf->GetPvsSize() < mTermMinPvs)
1070                ++ mStat.minPvsNodes;
1071
1072        mStat.accPvsSize += leaf->GetPvsSize();
1073
1074        if ((int)leaf->GetRays().size() < mTermMinRays)
1075                ++ mStat.minRaysNodes;
1076
1077        if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
1078                ++ mStat.maxRayContribNodes;
1079
1080        if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
1081                ++ mStat.minSizeNodes;
1082}
1083
1084
1085inline bool VspKdTree::TerminationCriteriaMet(const VspKdLeaf *leaf,
1086                                                                                          const AxisAlignedBox3 &box) const
1087{
1088        return ((leaf->GetPvsSize() < mTermMinPvs) ||
1089                    (leaf->mRays.size() < mTermMinRays) ||
1090                        (leaf->GetAvgRayContribution() > mTermMaxRayContribution ) ||
1091                        (leaf->mDepth >= mTermMaxDepth) ||
1092                        (SqrMagnitude(box.Size()) <= mTermMinSize));
1093}
1094
1095
1096VspKdNode *VspKdTree::SubdivideNode(VspKdLeaf *leaf,
1097                                                                        const AxisAlignedBox3 &box,
1098                                                                        AxisAlignedBox3 &backBBox,
1099                                                                        AxisAlignedBox3 &frontBBox)
1100{
1101        if (TerminationCriteriaMet(leaf, box))
1102        {
1103                if (1 && leaf->mDepth >= mTermMaxDepth)
1104                {
1105                        Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
1106                        Debug << "Bbox: " << GetBBox(leaf) << endl;
1107                }
1108                //Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
1109
1110                return leaf;
1111        }
1112
1113        float position;
1114        // first count ray sides
1115        int raysBack;
1116        int raysFront;
1117        int pvsBack;
1118        int pvsFront;
1119
1120        // select subdivision axis
1121        const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
1122        //Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" <<       pvsFront << endl;
1123
1124        if (axis == -1)
1125                return leaf;
1126       
1127        mStat.nodes += 2;
1128        ++ mStat.splits[axis];
1129
1130        // add the new nodes to the tree
1131        VspKdInterior *node = new VspKdInterior(leaf->mParent);
1132
1133        node->mAxis = axis;
1134        node->mPosition = position;
1135        node->mBox = box;
1136
1137        backBBox = box;
1138        frontBBox = box;
1139
1140        VspKdLeaf *back = new VspKdLeaf(node, raysBack, leaf->GetMaxCostMisses());
1141        back->SetPvsSize(pvsBack);
1142        VspKdLeaf *front = new VspKdLeaf(node, raysFront, leaf->GetMaxCostMisses());
1143        front->SetPvsSize(pvsFront);
1144
1145        // replace a link from node's parent
1146        if (leaf->mParent)
1147                leaf->mParent->ReplaceChildLink(leaf, node);
1148        // and setup child links
1149        node->SetupChildLinks(back, front);
1150
1151        if (axis <= VspKdNode::SPLIT_Z)
1152        {
1153                backBBox.SetMax(axis, position);
1154                frontBBox.SetMin(axis, position);
1155
1156                for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1157                                ri != leaf->mRays.end(); ++ ri)
1158                {
1159                        if ((*ri).mRay->IsActive())
1160                        {
1161                                // first unref ray from the former leaf
1162                                (*ri).mRay->Unref();
1163                                  float t;
1164                                // determine the side of this ray with respect to the plane
1165                                int side = node->ComputeRayIntersection(*ri, t);
1166
1167                                if (side == 0)
1168                                {
1169                                        if ((*ri).mRay->HasPosDir(axis))
1170                                        {
1171                                                back->AddRay(RayInfo((*ri).mRay,
1172                                                                                         (*ri).mMinT,
1173                                                                                         t));
1174                                                front->AddRay(RayInfo((*ri).mRay,
1175                                                                                          t,
1176                                                                                          (*ri).mMaxT));
1177                                        }
1178                                        else
1179                                        {
1180                                                back->AddRay(RayInfo((*ri).mRay,
1181                                                                                         t,
1182                                                                                         (*ri).mMaxT));
1183                                                front->AddRay(RayInfo((*ri).mRay,
1184                                                                                          (*ri).mMinT,
1185                                                                                          t));
1186                                        }
1187                                }
1188                                else
1189                                        if (side == 1)
1190                                                front->AddRay(*ri);
1191                                        else
1192                                                back->AddRay(*ri);
1193                        } else
1194                                (*ri).mRay->Unref();
1195                }
1196        }
1197        else
1198        {
1199                // rays front/back
1200
1201        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1202                        ri != leaf->mRays.end(); ++ ri)
1203                {
1204                        if ((*ri).mRay->IsActive())
1205                        {
1206                                // first unref ray from the former leaf
1207                                (*ri).mRay->Unref();
1208
1209                                int side;
1210                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1211                                        side = 1;
1212                                else
1213                                        side = -1;
1214
1215                                if (side == 1)
1216                                        front->AddRay(*ri);
1217                                else
1218                                        back->AddRay(*ri);
1219                        }
1220                        else
1221                                (*ri).mRay->Unref();
1222                }
1223        }
1224
1225        // update stats
1226        mStat.rayRefs -= (int)leaf->mRays.size();
1227        mStat.rayRefs += raysBack + raysFront;
1228
1229        DEL_PTR(leaf);
1230
1231        return node;
1232}
1233
1234int VspKdTree::ReleaseMemory(const int time)
1235{
1236        stack<VspKdNode *> tstack;
1237
1238        // find a node in the tree which subtree will be collapsed
1239        int maxAccessTime = time - mAccessTimeThreshold;
1240        int released = 0;
1241
1242        tstack.push(mRoot);
1243
1244        while (!tstack.empty())
1245        {
1246                VspKdNode *node = tstack.top();
1247                tstack.pop();
1248
1249                if (!node->IsLeaf())
1250                {
1251                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
1252                        // cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1253                        if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
1254                        {
1255                                released = CollapseSubtree(node, time);
1256                                break;
1257                        }
1258                        if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
1259                        {
1260                                tstack.push(in->GetFront());
1261                                tstack.push(in->GetBack());
1262                        }
1263                        else
1264                        {
1265                                tstack.push(in->GetBack());
1266                                tstack.push(in->GetFront());
1267                        }
1268                }
1269        }
1270
1271        while (tstack.empty())
1272        {
1273                // could find node to collaps...
1274                // cout<<"Could not find a node to release "<<endl;
1275                break;
1276        }
1277
1278        return released;
1279}
1280
1281
1282VspKdNode *VspKdTree::SubdivideLeaf(VspKdLeaf *leaf,
1283                                                                                const float sizeThreshold)
1284{
1285        VspKdNode *node = leaf;
1286
1287        AxisAlignedBox3 leafBBox = GetBBox(leaf);
1288
1289        static int pass = 0;
1290        ++ pass;
1291
1292        // check if we should perform a dynamic subdivision of the leaf
1293        if (// leaf->mRays.size() > (unsigned)termMinCost &&
1294                (leaf->GetPvsSize() >= mTermMinPvs) &&
1295                (SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
1296        {
1297        // memory check and release
1298                if (GetMemUsage() > mMaxTotalMemory)
1299                        ReleaseMemory(pass);
1300
1301                AxisAlignedBox3 backBBox, frontBBox;
1302
1303                // subdivide the node
1304                node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
1305        }
1306        return node;
1307}
1308
1309
1310
1311void VspKdTree::UpdateRays(VssRayContainer &remove,
1312                                                   VssRayContainer &add)
1313{
1314        VspKdLeaf::NewMail();
1315
1316        // schedule rays for removal
1317        for(VssRayContainer::const_iterator ri = remove.begin();
1318                ri != remove.end(); ++ ri)
1319        {
1320                (*ri)->ScheduleForRemoval();
1321        }
1322
1323        int inactive = 0;
1324
1325        for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
1326        {
1327                if ((*ri)->ScheduledForRemoval())
1328                        //      RemoveRay(*ri, NULL, false);
1329                        // !!! BUG - with true it does not work correctly - aggreated delete
1330                        RemoveRay(*ri, NULL, true);
1331                else
1332                        ++ inactive;
1333        }
1334
1335        //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1336        for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
1337        {
1338                AddRay(*ri);
1339        }
1340}
1341
1342
1343void VspKdTree::RemoveRay(VssRay *ray,
1344                                                  vector<VspKdLeaf *> *affectedLeaves,
1345                                                  const bool removeAllScheduledRays)
1346{
1347        stack<RayTraversalData> tstack;
1348
1349        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1350
1351        RayTraversalData data;
1352
1353        // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1354        while (!tstack.empty())
1355        {
1356                data = tstack.top();
1357                tstack.pop();
1358
1359                if (!data.mNode->IsLeaf())
1360                {
1361                        // split the set of rays in two groups intersecting the
1362                        // two subtrees
1363                        TraverseInternalNode(data, tstack);
1364        }
1365                else
1366                {
1367                        // remove the ray from the leaf
1368                        // find the ray in the leaf and swap it with the last ray...
1369                        VspKdLeaf *leaf = (VspKdLeaf *)data.mNode;
1370
1371                        if (!leaf->Mailed())
1372                        {
1373                                leaf->Mail();
1374
1375                                if (affectedLeaves)
1376                                        affectedLeaves->push_back(leaf);
1377
1378                                if (removeAllScheduledRays)
1379                                {
1380                                        int tail = (int)leaf->mRays.size() - 1;
1381
1382                                        for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
1383                                        {
1384                                                if (leaf->mRays[i].mRay->ScheduledForRemoval())
1385                                                {
1386                                                        // find a ray to replace it with
1387                                                        while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
1388                                                        {
1389                                                                ++ mStat.removedRayRefs;
1390                                                                leaf->mRays[tail].mRay->Unref();
1391                                                                leaf->mRays.pop_back();
1392
1393                                                                -- tail;
1394                                                        }
1395
1396                                                        if (tail < i)
1397                                                                break;
1398
1399                                                        ++ mStat.removedRayRefs;
1400
1401                                                        leaf->mRays[i].mRay->Unref();
1402                                                        leaf->mRays[i] = leaf->mRays[tail];
1403                                                        leaf->mRays.pop_back();
1404                                                        -- tail;
1405                                                }
1406                                        }
1407                                }
1408                        }
1409
1410                        if (!removeAllScheduledRays)
1411                                for (int i=0; i < (int)leaf->mRays.size(); i++)
1412                                {
1413                                        if (leaf->mRays[i].mRay == ray)
1414                                        {
1415                                                ++ mStat.removedRayRefs;
1416                                                ray->Unref();
1417                                                leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
1418                                                leaf->mRays.pop_back();
1419                                            // check this ray again
1420                                                break;
1421                                        }
1422                                }
1423                }
1424        }
1425
1426        if (ray->RefCount() != 0)
1427        {
1428                cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1429                exit(1);
1430        }
1431
1432}
1433
1434void VspKdTree::AddRay(VssRay *ray)
1435{
1436        stack<RayTraversalData> tstack;
1437
1438        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1439
1440        RayTraversalData data;
1441
1442        while (!tstack.empty())
1443        {
1444                data = tstack.top();
1445                tstack.pop();
1446
1447                if (!data.mNode->IsLeaf())
1448                {
1449                        TraverseInternalNode(data, tstack);
1450                }
1451                else
1452                {
1453                        // remove the ray from the leaf
1454                        // find the ray in the leaf and swap it with the last ray
1455                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(data.mNode);
1456
1457                        leaf->AddRay(data.mRayData);
1458                        ++ mStat.addedRayRefs;
1459                }
1460        }
1461}
1462
1463void VspKdTree::TraverseInternalNode(RayTraversalData &data,
1464                                                                         stack<RayTraversalData> &tstack)
1465{
1466        VspKdInterior *in = dynamic_cast<VspKdInterior *>(data.mNode);
1467
1468        if (in->mAxis <= VspKdNode::SPLIT_Z)
1469        {
1470            // determine the side of this ray with respect to the plane
1471                float t;
1472                const int side = in->ComputeRayIntersection(data.mRayData, t);
1473
1474                if (side == 0)
1475                {
1476                        if (data.mRayData.mRay->HasPosDir(in->mAxis))
1477                        {
1478                                tstack.push(RayTraversalData(in->GetBack(),
1479                                                        RayInfo(data.mRayData.mRay, data.mRayData.mMinT, t)));
1480
1481                                tstack.push(RayTraversalData(in->GetFront(),
1482                                                        RayInfo(data.mRayData.mRay, t, data.mRayData.mMaxT)));
1483
1484                        }
1485                        else
1486                        {
1487                                tstack.push(RayTraversalData(in->GetBack(),
1488                                                        RayInfo(data.mRayData.mRay,
1489                                                        t,
1490                                                        data.mRayData.mMaxT)));
1491                                tstack.push(RayTraversalData(in->GetFront(),
1492                                                        RayInfo(data.mRayData.mRay,
1493                                                        data.mRayData.mMinT,
1494                                                        t)));
1495                        }
1496                }
1497                else
1498                        if (side == 1)
1499                                tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1500                        else
1501                                tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1502        }
1503        else
1504        {
1505                // directional split
1506                if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
1507                        tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1508                else
1509                        tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1510        }
1511}
1512
1513
1514int VspKdTree::CollapseSubtree(VspKdNode *sroot, const int time)
1515{
1516        // first count all rays in the subtree
1517        // use mail 1 for this purpose
1518        stack<VspKdNode *> tstack;
1519
1520        int rayCount = 0;
1521        int totalRayCount = 0;
1522        int collapsedNodes = 0;
1523
1524#if DEBUG_COLLAPSE
1525        cout << "Collapsing subtree" << endl;
1526        cout << "accessTime=" << sroot->GetAccessTime() << endl;
1527        cout << "depth=" << (int)sroot->depth << endl;
1528#endif
1529
1530        // tstat.collapsedSubtrees++;
1531        // tstat.collapseDepths += (int)sroot->depth;
1532        // tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1533        tstack.push(sroot);
1534        VssRay::NewMail();
1535
1536        while (!tstack.empty())
1537        {
1538                ++ collapsedNodes;
1539
1540                VspKdNode *node = tstack.top();
1541                tstack.pop();
1542                if (node->IsLeaf())
1543                {
1544                        VspKdLeaf *leaf = (VspKdLeaf *) node;
1545
1546                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1547                                        ri != leaf->mRays.end(); ++ ri)
1548                        {
1549                                ++ totalRayCount;
1550
1551                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed())
1552                                {
1553                                        (*ri).mRay->Mail();
1554                                        ++ rayCount;
1555                                }
1556                        }
1557                }
1558                else
1559                {
1560                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetFront());
1561                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetBack());
1562                }
1563        }
1564
1565        VssRay::NewMail();
1566
1567        // create a new node that will hold the rays
1568        VspKdLeaf *newLeaf = new VspKdLeaf(sroot->mParent, rayCount);
1569
1570        if (newLeaf->mParent)
1571        {
1572                newLeaf->mParent->ReplaceChildLink(sroot, newLeaf);
1573        }
1574        tstack.push(sroot);
1575
1576        while (!tstack.empty())
1577        {
1578                VspKdNode *node = tstack.top();
1579                tstack.pop();
1580
1581                if (node->IsLeaf())
1582                {
1583                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
1584
1585                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1586                                        ri != leaf->mRays.end(); ++ ri)
1587                        {
1588                                // unref this ray from the old node
1589                                if ((*ri).mRay->IsActive())
1590                                {
1591                                        (*ri).mRay->Unref();
1592                                        if (!(*ri).mRay->Mailed())
1593                                        {
1594                                                (*ri).mRay->Mail();
1595                                                newLeaf->AddRay(*ri);
1596                                        }
1597                                }
1598                                else
1599                                        (*ri).mRay->Unref();
1600                        }
1601                }
1602                else
1603                {
1604                        VspKdInterior *interior =
1605                                dynamic_cast<VspKdInterior *>(node);
1606                        tstack.push(interior->GetBack());
1607                        tstack.push(interior->GetFront());
1608                }
1609        }
1610
1611        // delete the node and all its children
1612        DEL_PTR(sroot);
1613
1614        // for(VspKdNode::SRayContainer::iterator ri = newleaf->mRays.begin();
1615    //      ri != newleaf->mRays.end(); ++ ri)
1616        //     (*ri).ray->UnMail(2);
1617
1618#if DEBUG_COLLAPSE
1619        cout << "Total memory before=" << GetMemUsage() << endl;
1620#endif
1621
1622        mStat.nodes -= collapsedNodes - 1;
1623        mStat.rayRefs -= totalRayCount - rayCount;
1624
1625#if DEBUG_COLLAPSE
1626        cout << "collapsed nodes" << collapsedNodes << endl;
1627        cout << "collapsed rays" << totalRayCount - rayCount << endl;
1628        cout << "Total memory after=" << GetMemUsage() << endl;
1629        cout << "================================" << endl;
1630#endif
1631
1632        //  tstat.collapsedNodes += collapsedNodes;
1633        //  tstat.collapsedRays += totalRayCount - rayCount;
1634    return totalRayCount - rayCount;
1635}
1636
1637
1638int VspKdTree::GetPvsSize(VspKdNode *node, const AxisAlignedBox3 &box) const
1639{
1640        stack<VspKdNode *> tstack;
1641        tstack.push(mRoot);
1642
1643        Intersectable::NewMail();
1644        int pvsSize = 0;
1645
1646        while (!tstack.empty())
1647        {
1648                VspKdNode *node = tstack.top();
1649                tstack.pop();
1650
1651          if (node->IsLeaf())
1652                {
1653                        VspKdLeaf *leaf = (VspKdLeaf *)node;
1654                        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1655                                 ri != leaf->mRays.end(); ++ ri)
1656                        {
1657                                if ((*ri).mRay->IsActive())
1658                                {
1659                                        Intersectable *object;
1660#if BIDIRECTIONAL_RAY
1661                                        object = (*ri).mRay->mOriginObject;
1662
1663                                        if (object && !object->Mailed())
1664                                        {
1665                                                ++ pvsSize;
1666                                                object->Mail();
1667                                        }
1668#endif
1669                                        object = (*ri).mRay->mTerminationObject;
1670                                        if (object && !object->Mailed())
1671                                        {
1672                                                ++ pvsSize;
1673                                                object->Mail();
1674                                        }
1675                                }
1676                        }
1677                }
1678                else
1679                {
1680                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
1681
1682                        if (in->mAxis < 3)
1683                        {
1684                                if (box.Max(in->mAxis) >= in->mPosition)
1685                                        tstack.push(in->GetFront());
1686
1687                                if (box.Min(in->mAxis) <= in->mPosition)
1688                                        tstack.push(in->GetBack());
1689                        }
1690                        else
1691                        {
1692                                // both nodes for directional splits
1693                                tstack.push(in->GetFront());
1694                                tstack.push(in->GetBack());
1695                        }
1696                }
1697        }
1698
1699        return pvsSize;
1700}
1701
1702void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1703                                                                                         float &maxRayContribution,
1704                                                                                         float &avgRayContribution)
1705{
1706        stack<VspKdNode *> tstack;
1707        tstack.push(mRoot);
1708
1709        minRayContribution = 1.0f;
1710        maxRayContribution = 0.0f;
1711
1712        float sumRayContribution = 0.0f;
1713        int leaves = 0;
1714
1715        while (!tstack.empty())
1716        {
1717                VspKdNode *node = tstack.top();
1718                tstack.pop();
1719                if (node->IsLeaf())
1720                {
1721                        leaves++;
1722                        VspKdLeaf *leaf = (VspKdLeaf *)node;
1723                        float c = leaf->GetAvgRayContribution();
1724
1725                        if (c > maxRayContribution)
1726                                maxRayContribution = c;
1727                        if (c < minRayContribution)
1728                                minRayContribution = c;
1729                        sumRayContribution += c;
1730
1731                }
1732                else
1733                {
1734                        VspKdInterior *in = (VspKdInterior *)node;
1735                       
1736                        tstack.push(in->GetFront());
1737                        tstack.push(in->GetBack());
1738                }
1739        }
1740
1741        Debug << "sum=" << sumRayContribution << endl;
1742        Debug << "leaves=" << leaves << endl;
1743        avgRayContribution = sumRayContribution / (float)leaves;
1744}
1745
1746
1747int VspKdTree::GenerateRays(const float ratioPerLeaf,
1748                                                        SimpleRayContainer &rays)
1749{
1750        stack<VspKdNode *> tstack;
1751        tstack.push(mRoot);
1752
1753        while (!tstack.empty())
1754        {
1755                VspKdNode *node = tstack.top();
1756                tstack.pop();
1757
1758                if (node->IsLeaf())
1759                {
1760                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
1761
1762                        float c = leaf->GetAvgRayContribution();
1763                        int num = (int)(c*ratioPerLeaf + 0.5);
1764                        //                      cout<<num<<" ";
1765
1766                        for (int i=0; i < num; i++)
1767                        {
1768                                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1769                                /*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1770                                Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1771                                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1772                                rays.push_back(SimpleRay(origin, direction));*/
1773                        }
1774
1775                }
1776                else
1777                {
1778                        VspKdInterior *in =
1779                                dynamic_cast<VspKdInterior *>(node);
1780               
1781                        tstack.push(in->GetFront());
1782                        tstack.push(in->GetBack());
1783                }
1784        }
1785
1786        return (int)rays.size();
1787}
1788
1789
1790float VspKdTree::GetAvgPvsSize()
1791{
1792        stack<VspKdNode *> tstack;
1793        tstack.push(mRoot);
1794
1795        int sumPvs = 0;
1796        int leaves = 0;
1797
1798        while (!tstack.empty())
1799        {
1800                VspKdNode *node = tstack.top();
1801                tstack.pop();
1802
1803                if (node->IsLeaf())
1804                {
1805                        VspKdLeaf *leaf = (VspKdLeaf *)node;
1806                        // update pvs size
1807                        leaf->UpdatePvsSize();
1808                        sumPvs += leaf->GetPvsSize();
1809                        leaves++;
1810                }
1811                else
1812                {
1813                        VspKdInterior *in = (VspKdInterior *)node;
1814                       
1815                        tstack.push(in->GetFront());
1816                        tstack.push(in->GetBack());
1817                }
1818        }
1819
1820        return sumPvs / (float)leaves;
1821}
1822
1823VspKdNode *VspKdTree::GetRoot() const
1824{
1825        return mRoot;
1826}
1827
1828AxisAlignedBox3 VspKdTree::GetBBox(VspKdNode *node) const
1829{
1830        if (node->mParent == NULL)
1831                return mBox;
1832
1833        if (!node->IsLeaf())
1834                return (dynamic_cast<VspKdInterior *>(node))->mBox;
1835
1836        if (node->mParent->mAxis >= 3)
1837                return node->mParent->mBox;
1838
1839        AxisAlignedBox3 box(node->mParent->mBox);
1840        if (node->mParent->GetFront() == node)
1841                box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
1842        else
1843                box.SetMax(node->mParent->mAxis, node->mParent->mPosition);
1844
1845        return box;
1846}
1847
1848int     VspKdTree::GetRootPvsSize() const
1849{
1850        return GetPvsSize(mRoot, mBox);
1851}
1852
1853const VspKdStatistics &VspKdTree::GetStatistics() const
1854{
1855        return mStat;
1856}
1857
1858void VspKdTree::AddRays(VssRayContainer &add)
1859{
1860        VssRayContainer remove;
1861        UpdateRays(remove, add);
1862}
1863
1864// return memory usage in MB
1865float VspKdTree::GetMemUsage() const
1866{
1867        return
1868                (sizeof(VspKdTree) +
1869                 mStat.Leaves() * sizeof(VspKdLeaf) +
1870                 mStat.Interior() * sizeof(VspKdInterior) +
1871                 mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
1872}
1873
1874float VspKdTree::GetRayMemUsage() const
1875{
1876        return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
1877}
1878
1879
1880int VspKdTree::ComputePvsSize(VspKdNode *node,
1881                                                          const RayInfoContainer &globalRays) const
1882{
1883        int pvsSize = 0;
1884
1885        const AxisAlignedBox3 box = GetBBox(node);
1886
1887        RayInfoContainer::const_iterator it, it_end = globalRays.end();
1888
1889        // TODO
1890
1891        // warning: implicit conversion from VssRay to Ray
1892        for (it = globalRays.begin(); it != globalRays.end(); ++ it)
1893                pvsSize += box.GetMinMaxT(*(*it).mRay, NULL, NULL);
1894
1895        return pvsSize;
1896}
1897
1898
1899void VspKdTree::CollectLeaves(vector<VspKdLeaf *> &leaves) const
1900{
1901        stack<VspKdNode *> nodeStack;
1902        nodeStack.push(mRoot);
1903
1904        while (!nodeStack.empty())
1905        {
1906                VspKdNode *node = nodeStack.top();
1907
1908                nodeStack.pop();
1909
1910                if (node->IsLeaf())
1911                {
1912                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
1913                        leaves.push_back(leaf);
1914                }
1915                else
1916                {
1917                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
1918
1919                        nodeStack.push(interior->GetBack());
1920                        nodeStack.push(interior->GetFront());
1921                }
1922        }
1923}
1924
1925int VspKdTree::FindNeighbors(VspKdLeaf *n,
1926                                                         vector<VspKdLeaf *> &neighbors,
1927                                                         bool onlyUnmailed)
1928{
1929        stack<VspKdNode *> nodeStack;
1930
1931        nodeStack.push(mRoot);
1932
1933        AxisAlignedBox3 box = GetBBox(n);
1934
1935        while (!nodeStack.empty())
1936        {
1937                VspKdNode *node = nodeStack.top();
1938                nodeStack.pop();
1939
1940                if (node->IsLeaf())
1941                {
1942                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
1943
1944                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
1945                                neighbors.push_back(leaf);
1946                }
1947                else
1948                {
1949                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
1950
1951                        if (interior->mPosition > box.Max(interior->mAxis))
1952                                nodeStack.push(interior->mBack);
1953                        else
1954                        {
1955                                if (interior->mPosition < box.Min(interior->mAxis))
1956                                        nodeStack.push(interior->mFront);
1957                                else
1958                                {
1959                                        // random decision
1960                                        nodeStack.push(interior->mBack);
1961                                        nodeStack.push(interior->mFront);
1962                                }
1963                        }
1964                }
1965        }
1966
1967        return (int)neighbors.size();
1968}
1969
1970
1971void VspKdTree::SetViewCellsManager(ViewCellsManager *vcm)
1972{
1973        mViewCellsManager = vcm;
1974}
1975
1976
1977int VspKdTree::CastLineSegment(const Vector3 &origin,
1978                                                           const Vector3 &termination,
1979                                                           vector<ViewCell *> &viewcells)
1980{
1981        int hits = 0;
1982
1983        float mint = 0.0f, maxt = 1.0f;
1984        const Vector3 dir = termination - origin;
1985
1986        stack<LineTraversalData> tStack;
1987
1988        Intersectable::NewMail();
1989
1990        Vector3 entp = origin;
1991        Vector3 extp = termination;
1992
1993        VspKdNode *node = mRoot;
1994        VspKdNode *farChild;
1995
1996        float position;
1997        int axis;
1998
1999        while (1)
2000        {
2001                if (!node->IsLeaf())
2002                {
2003                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
2004                        position = in->mPosition;
2005                        axis = in->mAxis;
2006
2007                        if (entp[axis] <= position)
2008                        {
2009                                if (extp[axis] <= position)
2010                                {
2011                                        node = in->mBack;
2012                                        // cases N1,N2,N3,P5,Z2,Z3
2013                                        continue;
2014                                } else
2015                                {
2016                                        // case N4
2017                                        node = in->mBack;
2018                                        farChild = in->mFront;
2019                                }
2020                        }
2021                        else
2022                        {
2023                                if (position <= extp[axis])
2024                                {
2025                                        node = in->mFront;
2026                                        // cases P1,P2,P3,N5,Z1
2027                                        continue;
2028                                }
2029                                else
2030                                {
2031                                        node = in->mFront;
2032                                        farChild = in->mBack;
2033                                        // case P4
2034                                }
2035                        }
2036
2037                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2038                        // case N4 or P4
2039                        float tdist = (position - origin[axis]) / dir[axis];
2040                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2041                        extp = origin + dir * tdist;
2042                        maxt = tdist;
2043                }
2044                else
2045                {
2046                        // compute intersection with all objects in this leaf
2047                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2048                        ViewCell *vc = leaf->GetViewCell();
2049
2050                        if (!vc->Mailed())
2051                        {
2052                                vc->Mail();
2053                                viewcells.push_back(vc);
2054                                ++ hits;
2055                        }
2056
2057                        if (0) //matt TODO: REMOVE LATER
2058                                leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2059
2060                        // get the next node from the stack
2061                        if (tStack.empty())
2062                                break;
2063
2064                        entp = extp;
2065                        mint = maxt;
2066                       
2067                        LineTraversalData &s  = tStack.top();
2068                        node = s.mNode;
2069                        extp = s.mExitPoint;
2070                        maxt = s.mMaxT;
2071                        tStack.pop();
2072                }
2073        }
2074
2075        return hits;
2076}
2077
2078
2079void VspKdTree::GenerateViewCell(VspKdLeaf *leaf)
2080{
2081        RayInfoContainer::const_iterator it, it_end = leaf->GetRays().end();
2082
2083        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(mViewCellsManager->GenerateViewCell());
2084        leaf->SetViewCell(vc);
2085
2086        vc->SetVolume(GetBBox(leaf).GetVolume());
2087        vc->SetArea(GetBBox(leaf).SurfaceArea());
2088
2089        vc->mLeaves.push_back(leaf);
2090
2091        for (it = leaf->GetRays().begin(); it != it_end; ++ it)
2092        {
2093                VssRay *ray = (*it).mRay;
2094
2095                if (ray->mTerminationObject)
2096                  vc->GetPvs().AddSample(ray->mTerminationObject);
2097
2098                if (ray->mOriginObject)
2099                  vc->GetPvs().AddSample(ray->mOriginObject);
2100        }
2101}
2102
2103
2104bool VspKdTree::MergeViewCells(VspKdLeaf *l1, VspKdLeaf *l2)
2105{
2106        //-- change pointer to view cells of all leaves associated
2107        //-- with the previous view cells
2108        VspKdViewCell *fVc = l1->GetViewCell();
2109        VspKdViewCell *bVc = l2->GetViewCell();
2110
2111        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(
2112                mViewCellsManager->MergeViewCells(*fVc, *bVc));
2113
2114        // if merge was unsuccessful
2115        if (!vc) return false;
2116
2117        // set new size of view cell
2118        vc->SetVolume(fVc->GetVolume() + bVc->GetVolume());
2119        // new area
2120        vc->SetArea(fVc->GetArea() + bVc->GetArea());
2121
2122        vector<VspKdLeaf *> fLeaves = fVc->mLeaves;
2123        vector<VspKdLeaf *> bLeaves = bVc->mLeaves;
2124
2125        vector<VspKdLeaf *>::const_iterator it;
2126
2127        //-- change view cells of all the other leaves the view cell belongs to
2128        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
2129        {
2130                (*it)->SetViewCell(vc);
2131                vc->mLeaves.push_back(*it);
2132        }
2133
2134        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
2135        {
2136                (*it)->SetViewCell(vc);
2137                vc->mLeaves.push_back(*it);
2138        }
2139
2140        vc->Mail();
2141
2142        // clean up old view cells
2143        DEL_PTR(fVc);
2144        DEL_PTR(bVc);
2145
2146        return true;
2147}
2148
2149void VspKdTree::CollectMergeCandidates()
2150{
2151        MergeCandidate::sOverallCost = 0;
2152        vector<VspKdLeaf *> leaves;
2153
2154        // collect the leaves, e.g., the "voxels" that will build the view cells
2155        CollectLeaves(leaves);
2156
2157        VspKdLeaf::NewMail();
2158
2159        vector<VspKdLeaf *>::const_iterator it, it_end = leaves.end();
2160
2161        // generate view cells
2162        for (it = leaves.begin(); it != it_end; ++ it)
2163                GenerateViewCell(*it);
2164
2165        // find merge candidates and push them into queue
2166        for (it = leaves.begin(); it != it_end; ++ it)
2167        {
2168                VspKdLeaf *leaf = *it;
2169                ViewCell *vc = leaf->GetViewCell();
2170                // no leaf is part of two merge candidates
2171                leaf->Mail();
2172                MergeCandidate::sOverallCost +=
2173                        vc->GetVolume() * vc->GetPvs().GetSize();
2174                vector<VspKdLeaf *> neighbors;
2175                FindNeighbors(leaf, neighbors, true);
2176
2177                vector<VspKdLeaf *>::const_iterator nit,
2178                        nit_end = neighbors.end();
2179
2180                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
2181                {
2182                        mMergeQueue.push(MergeCandidate(leaf, *nit));
2183                }
2184        }
2185}
2186#if 0
2187void VspKdTree::CollectMergeCandidates(const vector<VspKdRay *> &rays)
2188{
2189        MergeCandidate::sOverallCost = 0;
2190
2191        vector<VspKdIntersection>::const_iterator iit;
2192        map<BspLeaf *, vector<BspLeaf*> candidateMap;
2193
2194        for (int i = 0; i < (int)rays.size(); ++ i)
2195        { 
2196                //VspKdLeaf::NewMail();
2197                VspKdRay *ray = rays[i];
2198         
2199                // traverse leaves stored in the rays and compare and merge consecutive
2200                // leaves (i.e., the neighbors in the tree)
2201                if (ray->intersections.size() < 2)
2202                        continue;
2203         
2204                iit = ray->intersections.begin();
2205
2206                BspLeaf *previousLeaf = (*iit).mLeaf;
2207       
2208                ++ iit;
2209               
2210                for (; iit != ray->intersections.end(); ++ iit)
2211                {
2212            BspLeaf *leaf = (*iit).mLeaf;
2213                        leaf->Mail
2214                        // TODO: how to sort out doubles?
2215                        MergeCandidate mc = MergeCandidate(leaf, previousLeaf);
2216                        mMergeQueue.push(mc);
2217
2218                        MergeCandidate::sOverallCost += mc.GetLeaf1Cost();
2219                        MergeCandidate::sOverallCost += mc.GetLeaf2Cost();
2220
2221                        previousLeaf = leaf;
2222        }
2223        }
2224}
2225#endif
2226
2227int VspKdTree::MergeViewCells(const VssRayContainer &rays)
2228{
2229        CollectMergeCandidates();
2230
2231        int merged = 0;
2232
2233        int vcSize = mStat.Leaves();
2234        // use priority queue to merge leaves
2235        while (!mMergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
2236                   (mMergeQueue.top().GetMergeCost() <
2237                    mMergeMaxCostRatio * MergeCandidate::sOverallCost))
2238        {
2239                //Debug << "mergecost: " << mergeQueue.top().GetMergeCost() /
2240                //MergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
2241                MergeCandidate mc = mMergeQueue.top();
2242                mMergeQueue.pop();
2243
2244                // both view cells equal!
2245                if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
2246                        continue;
2247
2248                if (mc.Valid())
2249                {
2250                        ViewCell::NewMail();
2251                        MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
2252
2253                        ++ merged;
2254                        -- vcSize;
2255                        // increase absolute merge cost
2256                        MergeCandidate::sOverallCost += mc.GetMergeCost();
2257                }
2258                // merge candidate not valid, because one of the leaves was already
2259                // merged with another one
2260                else
2261                {
2262                        // validate and reinsert into queue
2263                        mc.SetValid();
2264                        mMergeQueue.push(mc);
2265                        //Debug << "validate " << mc.GetMergeCost() << endl;
2266                }
2267        }
2268
2269        Debug << "merged " << merged << " of " << mStat.Leaves() << " leaves" << endl;
2270
2271        //TODO: should return sample contributions
2272        return merged;
2273}
2274
2275
2276void VspKdTree::RepairVcLeafLists()
2277{
2278        // list not valid anymore => clear
2279        stack<VspKdNode *> nodeStack;
2280        nodeStack.push(mRoot);
2281
2282        ViewCell::NewMail();
2283
2284        while (!nodeStack.empty())
2285        {
2286                VspKdNode *node = nodeStack.top();
2287                nodeStack.pop();
2288
2289                if (node->IsLeaf())
2290                {
2291                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2292
2293                        VspKdViewCell *viewCell = leaf->GetViewCell();
2294
2295                        if (!viewCell->Mailed())
2296                        {
2297                                viewCell->mLeaves.clear();
2298                                viewCell->Mail();
2299                        }
2300
2301                        viewCell->mLeaves.push_back(leaf);
2302                }
2303                else
2304                {
2305                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2306
2307                        nodeStack.push(interior->GetFront());
2308                        nodeStack.push(interior->GetBack());
2309                }
2310        }
2311}
2312
2313
2314VspKdNode *VspKdTree::CollapseTree(VspKdNode *node)
2315{
2316    if (node->IsLeaf())
2317                return node;
2318
2319        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2320
2321        VspKdNode *front = CollapseTree(interior->GetFront());
2322        VspKdNode *back = CollapseTree(interior->GetBack());
2323
2324        if (front->IsLeaf() && back->IsLeaf())
2325        {
2326                VspKdLeaf *frontLeaf = dynamic_cast<VspKdLeaf *>(front);
2327                VspKdLeaf *backLeaf = dynamic_cast<VspKdLeaf *>(back);
2328
2329                //-- collapse tree
2330                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
2331                {
2332                        VspKdViewCell *vc = frontLeaf->GetViewCell();
2333
2334                        VspKdLeaf *leaf = new VspKdLeaf(interior->GetParent(), 0);
2335                        leaf->SetViewCell(vc);
2336
2337                        // replace a link from node's parent
2338                        if (leaf->mParent)
2339                                leaf->mParent->ReplaceChildLink(node, leaf);
2340
2341                        delete interior;
2342
2343                        return leaf;
2344                }
2345        }
2346
2347        // revalidate leaves
2348        RepairVcLeafLists();
2349
2350        return node;
2351}
2352
2353
2354int VspKdTree::RefineViewCells(const VssRayContainer &rays)
2355{
2356        int shuffled = 0;
2357
2358        Debug << "refining " << (int)mMergeQueue.size() << " candidates " << endl;
2359        BspLeaf::NewMail();
2360
2361        // Use priority queue of remaining leaf pairs
2362        // These candidates either share the same view cells or
2363        // are border leaves which share a boundary.
2364        // We test if they can be shuffled, i.e.,
2365        // either one leaf is made part of one view cell or the other
2366        // leaf is made part of the other view cell. It is tested if the
2367        // remaining view cells are "better" than the old ones.
2368        while (!mMergeQueue.empty())
2369        {
2370                MergeCandidate mc = mMergeQueue.top();
2371                mMergeQueue.pop();
2372
2373                // both view cells equal or already shuffled
2374                if ((mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell()) ||
2375                        (mc.GetLeaf1()->Mailed()) || (mc.GetLeaf2()->Mailed()))
2376                        continue;
2377               
2378                // candidate for shuffling
2379                const bool wasShuffled = false;
2380                        //ShuffleLeaves(mc.GetLeaf1(), mc.GetLeaf2());
2381               
2382                //-- stats
2383                if (wasShuffled)
2384                        ++ shuffled;
2385        }
2386
2387        return shuffled;
2388}
2389
2390
2391inline int AddedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2392{
2393        return pvs1.AddPvs(pvs2);
2394}
2395
2396
2397inline int SubtractedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2398{
2399        return pvs1.SubtractPvs(pvs2);
2400}
2401
2402
2403float GetShuffledVcCost(VspKdLeaf *leaf, VspKdViewCell *vc1, VspKdViewCell *vc2)
2404{
2405#if 0 // TODO
2406        const int pvs1 = SubtractedPvsSize(vc1->GetPvs(), *leaf->mPvs);
2407        const int pvs2 = AddedPvsSize(vc2->GetPvs(), *leaf->mPvs);
2408
2409        const float area1 = vc1->GetArea() - leaf->mArea;
2410        const float area2 = vc2->GetArea() + leaf->mArea;
2411
2412        const float cost1 = pvs1 * area1;
2413        const float cost2 = pvs2 * area2;
2414
2415        return cost1 + cost2;
2416#else
2417        return -1;
2418#endif
2419}
2420
2421
2422void VspKdTree::ShuffleLeaf(VspKdLeaf *leaf,
2423                                                        VspKdViewCell *vc1,
2424                                                        VspKdViewCell *vc2) const
2425{
2426#if 0 // TODO
2427        // compute new pvs and area
2428        vc1->GetPvs().SubtractPvs(*leaf->mPvs);
2429        vc2->GetPvs().AddPvs(*leaf->mPvs);
2430       
2431        vc1->SetArea(vc1->GetArea() - leaf->mArea);
2432        vc2->SetArea(vc2->GetArea() + leaf->mArea);
2433
2434        /// add to second view cell
2435        vc2->mLeaves.push_back(leaf);
2436
2437        // erase leaf from old view cell
2438        vector<VspKdLeaf *>::iterator it = vc1->mLeaves.begin();
2439
2440        for (; *it != leaf; ++ it);
2441        vc1->mLeaves.erase(it);
2442
2443        leaf->SetViewCell(vc2);  // finally change view cell
2444#endif
2445}
2446
2447
2448bool VspKdTree::ShuffleLeaves(VspKdLeaf *leaf1, VspKdLeaf *leaf2) const
2449{
2450        VspKdViewCell *vc1 = leaf1->GetViewCell();
2451        VspKdViewCell *vc2 = leaf2->GetViewCell();
2452
2453        const float cost1 = vc1->GetPvs().GetSize() * vc1->GetArea();
2454        const float cost2 = vc2->GetPvs().GetSize() * vc2->GetArea();
2455
2456        const float oldCost = cost1 + cost2;
2457       
2458        float shuffledCost1 = Limits::Infinity;
2459        float shuffledCost2 = Limits::Infinity;
2460
2461        // the view cell should not be empty after the shuffle
2462        if (vc1->mLeaves.size() > 1)
2463                shuffledCost1 = GetShuffledVcCost(leaf1, vc1, vc2);
2464        if (vc2->mLeaves.size() > 1)
2465                shuffledCost2 = GetShuffledVcCost(leaf2, vc2, vc1);
2466
2467        // shuffling unsuccessful
2468        if ((oldCost <= shuffledCost1) && (oldCost <= shuffledCost2))
2469                return false;
2470       
2471        if (shuffledCost1 < shuffledCost2)
2472        {
2473                ShuffleLeaf(leaf1, vc1, vc2);
2474                leaf1->Mail();
2475        }
2476        else
2477        {
2478                ShuffleLeaf(leaf2, vc2, vc1);
2479                leaf2->Mail();
2480        }
2481
2482        return true;
2483}
2484
2485
2486
2487
2488/*********************************************************************/
2489/*                MergeCandidate implementation                      */
2490/*********************************************************************/
2491
2492
2493MergeCandidate::MergeCandidate(VspKdLeaf *l1, VspKdLeaf *l2):
2494mMergeCost(0),
2495mLeaf1(l1),
2496mLeaf2(l2),
2497mLeaf1Id(l1->GetViewCell()->mMailbox),
2498mLeaf2Id(l2->GetViewCell()->mMailbox)
2499{
2500        EvalMergeCost();
2501}
2502
2503float MergeCandidate::GetLeaf1Cost() const
2504{
2505        VspKdViewCell *vc = mLeaf1->GetViewCell();
2506        return vc->GetPvs().GetSize() * vc->GetVolume();
2507}
2508
2509float MergeCandidate::GetLeaf2Cost() const
2510{
2511        VspKdViewCell *vc = mLeaf2->GetViewCell();
2512        return vc->GetPvs().GetSize() * vc->GetVolume();
2513}
2514
2515void MergeCandidate::EvalMergeCost()
2516{
2517        //-- compute pvs difference
2518        VspKdViewCell *vc1 = mLeaf1->GetViewCell();
2519        VspKdViewCell *vc2 = mLeaf2->GetViewCell();
2520
2521        const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
2522        const int vcPvs = diff1 + vc1->GetPvs().GetSize();
2523
2524        //-- compute ratio of old cost
2525        //-- (added size of left and right view cell times pvs size)
2526        //-- to new rendering cost (size of merged view cell times pvs size)
2527        const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();
2528       
2529        const float newCost =
2530                (float)vcPvs * (vc1->GetVolume() + vc2->GetVolume());
2531
2532        mMergeCost = newCost - oldCost;
2533
2534//      if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
2535//              mMergeCost += 1.0;
2536}
2537
2538void MergeCandidate::SetLeaf1(VspKdLeaf *l)
2539{
2540        mLeaf1 = l;
2541}
2542
2543void MergeCandidate::SetLeaf2(VspKdLeaf *l)
2544{
2545        mLeaf2 = l;
2546}
2547
2548
2549VspKdLeaf *MergeCandidate::GetLeaf1()
2550{
2551        return mLeaf1;
2552}
2553
2554
2555VspKdLeaf *MergeCandidate::GetLeaf2()
2556{
2557        return mLeaf2;
2558}
2559
2560
2561bool MergeCandidate::Valid() const
2562{
2563        return
2564                (mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
2565                (mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
2566}
2567
2568
2569float MergeCandidate::GetMergeCost() const
2570{
2571        return mMergeCost;
2572}
2573
2574
2575void MergeCandidate::SetValid()
2576{
2577        mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
2578        mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;
2579
2580        EvalMergeCost();
2581}
Note: See TracBrowser for help on using the repository browser.