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

Revision 582, 60.3 KB checked in by mattausch, 19 years ago (diff)

fixed bug in mergueue to find root of merge and sort out doube view cells

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