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

Revision 667, 59.5 KB checked in by mattausch, 19 years ago (diff)

test version
build script for testing. code is set uo for testing also

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