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

Revision 744, 59.5 KB checked in by mattausch, 18 years ago (diff)
Line 
1
2// ================================================================
3// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4// ****************************************************************
5/**
6   The KD tree based LSDS
7*/
8// Initial coding by
9/**
10   @author Jiri Bittner
11*/
12
13// Standard headers
14#include <stack>
15#include <queue>
16#include <algorithm>
17#include <fstream>
18#include <string>
19
20#include "VspKdTree.h"
21
22#include "Environment.h"
23#include "VssRay.h"
24#include "Intersectable.h"
25#include "Ray.h"
26#include "RayInfo.h"
27#include "ViewCellsManager.h"
28#include "ViewCellBsp.h"
29
30// Static variables
31int VspKdLeaf::sMailId = 0;
32int 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        {
1128                ++ mStat.maxDepthNodes;
1129        }
1130
1131        if (leaf->GetPvsSize() < mTermMinPvs)
1132                ++ mStat.minPvsNodes;
1133
1134        mStat.accPvsSize += leaf->GetPvsSize();
1135
1136        if ((int)leaf->GetRays().size() < mTermMinRays)
1137                ++ mStat.minRaysNodes;
1138
1139        if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
1140                ++ mStat.maxRayContribNodes;
1141
1142        if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
1143                ++ mStat.minSizeNodes;
1144}
1145
1146
1147inline bool VspKdTree::TerminationCriteriaMet(const VspKdLeaf *leaf,
1148                                                                                          const AxisAlignedBox3 &box) const
1149{
1150        return ((leaf->GetPvsSize() < mTermMinPvs) ||
1151                    (leaf->mRays.size() < mTermMinRays) ||
1152                        (leaf->GetAvgRayContribution() > mTermMaxRayContribution ) ||
1153                        (leaf->mDepth >= mTermMaxDepth) ||
1154                        (mStat.Leaves() >= mMaxViewCells) ||
1155                        (SqrMagnitude(box.Size()) <= mTermMinSize));
1156}
1157
1158
1159VspKdNode *VspKdTree::SubdivideNode(VspKdLeaf *leaf,
1160                                                                        const AxisAlignedBox3 &box,
1161                                                                        AxisAlignedBox3 &backBBox,
1162                                                                        AxisAlignedBox3 &frontBBox)
1163{
1164        if (TerminationCriteriaMet(leaf, box))
1165        {
1166                if (1 && leaf->mDepth >= mTermMaxDepth)
1167                {
1168                        Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
1169                        Debug << "Bbox: " << GetBBox(leaf) << endl;
1170                }
1171                //Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
1172
1173                return leaf;
1174        }
1175
1176        float position;
1177        // first count ray sides
1178        int raysBack;
1179        int raysFront;
1180        int pvsBack;
1181        int pvsFront;
1182
1183        // select subdivision axis
1184        const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
1185        //Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" <<       pvsFront << endl;
1186
1187        if (axis == -1)
1188                return leaf;
1189       
1190        mStat.nodes += 2;
1191        ++ mStat.splits[axis];
1192
1193        // add the new nodes to the tree
1194        VspKdInterior *node =
1195                new VspKdInterior(dynamic_cast<VspKdInterior *>(leaf->mParent));
1196
1197        node->mAxis = axis;
1198        node->mPosition = position;
1199        node->mBox = box;
1200
1201        backBBox = box;
1202        frontBBox = box;
1203
1204        VspKdLeaf *back = new VspKdLeaf(node, raysBack, leaf->GetMaxCostMisses());
1205        back->SetPvsSize(pvsBack);
1206        VspKdLeaf *front = new VspKdLeaf(node, raysFront, leaf->GetMaxCostMisses());
1207        front->SetPvsSize(pvsFront);
1208
1209        // replace a link from node's parent
1210        if (leaf->mParent)
1211        {
1212                VspKdInterior *parent = dynamic_cast<VspKdInterior *>(leaf->mParent);
1213                parent->ReplaceChildLink(leaf, node);
1214        }
1215
1216        // and setup child links
1217        node->SetupChildLinks(back, front);
1218
1219        if (axis <= VspKdNode::SPLIT_Z)
1220        {
1221                backBBox.SetMax(axis, position);
1222                frontBBox.SetMin(axis, position);
1223
1224                for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1225                                ri != leaf->mRays.end(); ++ ri)
1226                {
1227                        if ((*ri).mRay->IsActive())
1228                        {
1229                                // first unref ray from the former leaf
1230                                (*ri).mRay->Unref();
1231                                  float t;
1232                                // determine the side of this ray with respect to the plane
1233                                int side = node->ComputeRayIntersection(*ri, t);
1234
1235                                if (side == 0)
1236                                {
1237                                        if ((*ri).mRay->HasPosDir(axis))
1238                                        {
1239                                                back->AddRay(RayInfo((*ri).mRay,
1240                                                                                         (*ri).mMinT,
1241                                                                                         t));
1242                                                front->AddRay(RayInfo((*ri).mRay,
1243                                                                                          t,
1244                                                                                          (*ri).mMaxT));
1245                                        }
1246                                        else
1247                                        {
1248                                                back->AddRay(RayInfo((*ri).mRay,
1249                                                                                         t,
1250                                                                                         (*ri).mMaxT));
1251                                                front->AddRay(RayInfo((*ri).mRay,
1252                                                                                          (*ri).mMinT,
1253                                                                                          t));
1254                                        }
1255                                }
1256                                else
1257                                        if (side == 1)
1258                                                front->AddRay(*ri);
1259                                        else
1260                                                back->AddRay(*ri);
1261                        } else
1262                                (*ri).mRay->Unref();
1263                }
1264        }
1265        else
1266        {
1267                // rays front/back
1268
1269        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1270                        ri != leaf->mRays.end(); ++ ri)
1271                {
1272                        if ((*ri).mRay->IsActive())
1273                        {
1274                                // first unref ray from the former leaf
1275                                (*ri).mRay->Unref();
1276
1277                                int side;
1278                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1279                                        side = 1;
1280                                else
1281                                        side = -1;
1282
1283                                if (side == 1)
1284                                        front->AddRay(*ri);
1285                                else
1286                                        back->AddRay(*ri);
1287                        }
1288                        else
1289                                (*ri).mRay->Unref();
1290                }
1291        }
1292
1293        // update stats
1294        mStat.rayRefs -= (int)leaf->mRays.size();
1295        mStat.rayRefs += raysBack + raysFront;
1296
1297        DEL_PTR(leaf);
1298
1299        return node;
1300}
1301
1302int VspKdTree::ReleaseMemory(const int time)
1303{
1304        stack<VspKdNode *> tstack;
1305
1306        // find a node in the tree which subtree will be collapsed
1307        int maxAccessTime = time - mAccessTimeThreshold;
1308        int released = 0;
1309
1310        tstack.push(mRoot);
1311
1312        while (!tstack.empty())
1313        {
1314                VspKdNode *node = tstack.top();
1315                tstack.pop();
1316
1317                if (!node->IsLeaf())
1318                {
1319                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
1320                        // cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1321                        if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
1322                        {
1323                                released = CollapseSubtree(node, time);
1324                                break;
1325                        }
1326                        if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
1327                        {
1328                                tstack.push(in->GetFront());
1329                                tstack.push(in->GetBack());
1330                        }
1331                        else
1332                        {
1333                                tstack.push(in->GetBack());
1334                                tstack.push(in->GetFront());
1335                        }
1336                }
1337        }
1338
1339        while (tstack.empty())
1340        {
1341                // could find node to collaps...
1342                // cout<<"Could not find a node to release "<<endl;
1343                break;
1344        }
1345
1346        return released;
1347}
1348
1349
1350VspKdNode *VspKdTree::SubdivideLeaf(VspKdLeaf *leaf,
1351                                                                                const float sizeThreshold)
1352{
1353        VspKdNode *node = leaf;
1354
1355        AxisAlignedBox3 leafBBox = GetBBox(leaf);
1356
1357        static int pass = 0;
1358        ++ pass;
1359
1360        // check if we should perform a dynamic subdivision of the leaf
1361        if (// leaf->mRays.size() > (unsigned)termMinCost &&
1362                (leaf->GetPvsSize() >= mTermMinPvs) &&
1363                (SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
1364        {
1365        // memory check and release
1366                if (GetMemUsage() > mMaxTotalMemory)
1367                        ReleaseMemory(pass);
1368
1369                AxisAlignedBox3 backBBox, frontBBox;
1370
1371                // subdivide the node
1372                node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
1373        }
1374        return node;
1375}
1376
1377
1378
1379void VspKdTree::UpdateRays(VssRayContainer &remove,
1380                                                   VssRayContainer &add)
1381{
1382        VspKdLeaf::NewMail();
1383
1384        // schedule rays for removal
1385        for(VssRayContainer::const_iterator ri = remove.begin();
1386                ri != remove.end(); ++ ri)
1387        {
1388                (*ri)->ScheduleForRemoval();
1389        }
1390
1391        int inactive = 0;
1392
1393        for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
1394        {
1395                if ((*ri)->ScheduledForRemoval())
1396                        //      RemoveRay(*ri, NULL, false);
1397                        // !!! BUG - with true it does not work correctly - aggreated delete
1398                        RemoveRay(*ri, NULL, true);
1399                else
1400                        ++ inactive;
1401        }
1402
1403        //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1404        for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
1405        {
1406                AddRay(*ri);
1407        }
1408}
1409
1410
1411void VspKdTree::RemoveRay(VssRay *ray,
1412                                                  vector<VspKdLeaf *> *affectedLeaves,
1413                                                  const bool removeAllScheduledRays)
1414{
1415        stack<RayTraversalData> tstack;
1416
1417        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1418
1419        RayTraversalData data;
1420
1421        // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1422        while (!tstack.empty())
1423        {
1424                data = tstack.top();
1425                tstack.pop();
1426
1427                if (!data.mNode->IsLeaf())
1428                {
1429                        // split the set of rays in two groups intersecting the
1430                        // two subtrees
1431                        TraverseInternalNode(data, tstack);
1432        }
1433                else
1434                {
1435                        // remove the ray from the leaf
1436                        // find the ray in the leaf and swap it with the last ray...
1437                        VspKdLeaf *leaf = (VspKdLeaf *)data.mNode;
1438
1439                        if (!leaf->Mailed())
1440                        {
1441                                leaf->Mail();
1442
1443                                if (affectedLeaves)
1444                                        affectedLeaves->push_back(leaf);
1445
1446                                if (removeAllScheduledRays)
1447                                {
1448                                        int tail = (int)leaf->mRays.size() - 1;
1449
1450                                        for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
1451                                        {
1452                                                if (leaf->mRays[i].mRay->ScheduledForRemoval())
1453                                                {
1454                                                        // find a ray to replace it with
1455                                                        while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
1456                                                        {
1457                                                                ++ mStat.removedRayRefs;
1458                                                                leaf->mRays[tail].mRay->Unref();
1459                                                                leaf->mRays.pop_back();
1460
1461                                                                -- tail;
1462                                                        }
1463
1464                                                        if (tail < i)
1465                                                                break;
1466
1467                                                        ++ mStat.removedRayRefs;
1468
1469                                                        leaf->mRays[i].mRay->Unref();
1470                                                        leaf->mRays[i] = leaf->mRays[tail];
1471                                                        leaf->mRays.pop_back();
1472                                                        -- tail;
1473                                                }
1474                                        }
1475                                }
1476                        }
1477
1478                        if (!removeAllScheduledRays)
1479                                for (int i=0; i < (int)leaf->mRays.size(); i++)
1480                                {
1481                                        if (leaf->mRays[i].mRay == ray)
1482                                        {
1483                                                ++ mStat.removedRayRefs;
1484                                                ray->Unref();
1485                                                leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
1486                                                leaf->mRays.pop_back();
1487                                            // check this ray again
1488                                                break;
1489                                        }
1490                                }
1491                }
1492        }
1493
1494        if (ray->RefCount() != 0)
1495        {
1496                cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1497                exit(1);
1498        }
1499
1500}
1501
1502void VspKdTree::AddRay(VssRay *ray)
1503{
1504        stack<RayTraversalData> tstack;
1505
1506        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1507
1508        RayTraversalData data;
1509
1510        while (!tstack.empty())
1511        {
1512                data = tstack.top();
1513                tstack.pop();
1514
1515                if (!data.mNode->IsLeaf())
1516                {
1517                        TraverseInternalNode(data, tstack);
1518                }
1519                else
1520                {
1521                        // remove the ray from the leaf
1522                        // find the ray in the leaf and swap it with the last ray
1523                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(data.mNode);
1524
1525                        leaf->AddRay(data.mRayData);
1526                        ++ mStat.addedRayRefs;
1527                }
1528        }
1529}
1530
1531void VspKdTree::TraverseInternalNode(RayTraversalData &data,
1532                                                                         stack<RayTraversalData> &tstack)
1533{
1534        VspKdInterior *in = dynamic_cast<VspKdInterior *>(data.mNode);
1535
1536        if (in->mAxis <= VspKdNode::SPLIT_Z)
1537        {
1538            // determine the side of this ray with respect to the plane
1539                float t;
1540                const int side = in->ComputeRayIntersection(data.mRayData, t);
1541
1542                if (side == 0)
1543                {
1544                        if (data.mRayData.mRay->HasPosDir(in->mAxis))
1545                        {
1546                                tstack.push(RayTraversalData(in->GetBack(),
1547                                                        RayInfo(data.mRayData.mRay, data.mRayData.mMinT, t)));
1548
1549                                tstack.push(RayTraversalData(in->GetFront(),
1550                                                        RayInfo(data.mRayData.mRay, t, data.mRayData.mMaxT)));
1551
1552                        }
1553                        else
1554                        {
1555                                tstack.push(RayTraversalData(in->GetBack(),
1556                                                        RayInfo(data.mRayData.mRay,
1557                                                        t,
1558                                                        data.mRayData.mMaxT)));
1559                                tstack.push(RayTraversalData(in->GetFront(),
1560                                                        RayInfo(data.mRayData.mRay,
1561                                                        data.mRayData.mMinT,
1562                                                        t)));
1563                        }
1564                }
1565                else
1566                        if (side == 1)
1567                                tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1568                        else
1569                                tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1570        }
1571        else
1572        {
1573                // directional split
1574                if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
1575                        tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1576                else
1577                        tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1578        }
1579}
1580
1581
1582int VspKdTree::CollapseSubtree(VspKdNode *sroot, const int time)
1583{
1584        // first count all rays in the subtree
1585        // use mail 1 for this purpose
1586        stack<VspKdNode *> tstack;
1587
1588        int rayCount = 0;
1589        int totalRayCount = 0;
1590        int collapsedNodes = 0;
1591
1592#if DEBUG_COLLAPSE
1593        cout << "Collapsing subtree" << endl;
1594        cout << "accessTime=" << sroot->GetAccessTime() << endl;
1595        cout << "depth=" << (int)sroot->depth << endl;
1596#endif
1597
1598        // tstat.collapsedSubtrees++;
1599        // tstat.collapseDepths += (int)sroot->depth;
1600        // tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1601        tstack.push(sroot);
1602        VssRay::NewMail();
1603
1604        while (!tstack.empty())
1605        {
1606                ++ collapsedNodes;
1607
1608                VspKdNode *node = tstack.top();
1609                tstack.pop();
1610                if (node->IsLeaf())
1611                {
1612                        VspKdLeaf *leaf = (VspKdLeaf *) node;
1613
1614                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1615                                        ri != leaf->mRays.end(); ++ ri)
1616                        {
1617                                ++ totalRayCount;
1618
1619                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed())
1620                                {
1621                                        (*ri).mRay->Mail();
1622                                        ++ rayCount;
1623                                }
1624                        }
1625                }
1626                else
1627                {
1628                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetFront());
1629                        tstack.push(dynamic_cast<VspKdInterior *>(node)->GetBack());
1630                }
1631        }
1632
1633        VssRay::NewMail();
1634
1635        // create a new node that will hold the rays
1636        VspKdLeaf *newLeaf = new VspKdLeaf(sroot->mParent, rayCount);
1637
1638        VspKdInterior *parent =
1639                dynamic_cast<VspKdInterior *>(newLeaf->mParent);
1640
1641        if (parent)
1642        {
1643                parent->ReplaceChildLink(sroot, newLeaf);
1644        }
1645        tstack.push(sroot);
1646
1647        while (!tstack.empty())
1648        {
1649                VspKdNode *node = tstack.top();
1650                tstack.pop();
1651
1652                if (node->IsLeaf())
1653                {
1654                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
1655
1656                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1657                                        ri != leaf->mRays.end(); ++ ri)
1658                        {
1659                                // unref this ray from the old node
1660                                if ((*ri).mRay->IsActive())
1661                                {
1662                                        (*ri).mRay->Unref();
1663                                        if (!(*ri).mRay->Mailed())
1664                                        {
1665                                                (*ri).mRay->Mail();
1666                                                newLeaf->AddRay(*ri);
1667                                        }
1668                                }
1669                                else
1670                                        (*ri).mRay->Unref();
1671                        }
1672                }
1673                else
1674                {
1675                        VspKdInterior *interior =
1676                                dynamic_cast<VspKdInterior *>(node);
1677                        tstack.push(interior->GetBack());
1678                        tstack.push(interior->GetFront());
1679                }
1680        }
1681
1682        // delete the node and all its children
1683        DEL_PTR(sroot);
1684
1685        // for(VspKdNode::SRayContainer::iterator ri = newleaf->mRays.begin();
1686    //      ri != newleaf->mRays.end(); ++ ri)
1687        //     (*ri).ray->UnMail(2);
1688
1689#if DEBUG_COLLAPSE
1690        cout << "Total memory before=" << GetMemUsage() << endl;
1691#endif
1692
1693        mStat.nodes -= collapsedNodes - 1;
1694        mStat.rayRefs -= totalRayCount - rayCount;
1695
1696#if DEBUG_COLLAPSE
1697        cout << "collapsed nodes" << collapsedNodes << endl;
1698        cout << "collapsed rays" << totalRayCount - rayCount << endl;
1699        cout << "Total memory after=" << GetMemUsage() << endl;
1700        cout << "================================" << endl;
1701#endif
1702
1703        //  tstat.collapsedNodes += collapsedNodes;
1704        //  tstat.collapsedRays += totalRayCount - rayCount;
1705    return totalRayCount - rayCount;
1706}
1707
1708
1709int VspKdTree::GetPvsSize(VspKdNode *node, const AxisAlignedBox3 &box) const
1710{
1711        stack<VspKdNode *> tstack;
1712        tstack.push(mRoot);
1713
1714        Intersectable::NewMail();
1715        int pvsSize = 0;
1716
1717        while (!tstack.empty())
1718        {
1719                VspKdNode *node = tstack.top();
1720                tstack.pop();
1721
1722          if (node->IsLeaf())
1723                {
1724                        VspKdLeaf *leaf = (VspKdLeaf *)node;
1725                        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1726                                 ri != leaf->mRays.end(); ++ ri)
1727                        {
1728                                if ((*ri).mRay->IsActive())
1729                                {
1730                                        Intersectable *object;
1731#if BIDIRECTIONAL_RAY
1732                                        object = (*ri).mRay->mOriginObject;
1733
1734                                        if (object && !object->Mailed())
1735                                        {
1736                                                ++ pvsSize;
1737                                                object->Mail();
1738                                        }
1739#endif
1740                                        object = (*ri).mRay->mTerminationObject;
1741                                        if (object && !object->Mailed())
1742                                        {
1743                                                ++ pvsSize;
1744                                                object->Mail();
1745                                        }
1746                                }
1747                        }
1748                }
1749                else
1750                {
1751                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
1752
1753                        if (in->mAxis < 3)
1754                        {
1755                                if (box.Max(in->mAxis) >= in->mPosition)
1756                                        tstack.push(in->GetFront());
1757
1758                                if (box.Min(in->mAxis) <= in->mPosition)
1759                                        tstack.push(in->GetBack());
1760                        }
1761                        else
1762                        {
1763                                // both nodes for directional splits
1764                                tstack.push(in->GetFront());
1765                                tstack.push(in->GetBack());
1766                        }
1767                }
1768        }
1769
1770        return pvsSize;
1771}
1772
1773void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1774                                                                                         float &maxRayContribution,
1775                                                                                         float &avgRayContribution)
1776{
1777        stack<VspKdNode *> tstack;
1778        tstack.push(mRoot);
1779
1780        minRayContribution = 1.0f;
1781        maxRayContribution = 0.0f;
1782
1783        float sumRayContribution = 0.0f;
1784        int leaves = 0;
1785
1786        while (!tstack.empty())
1787        {
1788                VspKdNode *node = tstack.top();
1789                tstack.pop();
1790                if (node->IsLeaf())
1791                {
1792                        leaves++;
1793                        VspKdLeaf *leaf = (VspKdLeaf *)node;
1794                        float c = leaf->GetAvgRayContribution();
1795
1796                        if (c > maxRayContribution)
1797                                maxRayContribution = c;
1798                        if (c < minRayContribution)
1799                                minRayContribution = c;
1800                        sumRayContribution += c;
1801
1802                }
1803                else
1804                {
1805                        VspKdInterior *in = (VspKdInterior *)node;
1806                       
1807                        tstack.push(in->GetFront());
1808                        tstack.push(in->GetBack());
1809                }
1810        }
1811
1812        Debug << "sum=" << sumRayContribution << endl;
1813        Debug << "leaves=" << leaves << endl;
1814        avgRayContribution = sumRayContribution / (float)leaves;
1815}
1816
1817
1818int VspKdTree::GenerateRays(const float ratioPerLeaf,
1819                                                        SimpleRayContainer &rays)
1820{
1821        stack<VspKdNode *> tstack;
1822        tstack.push(mRoot);
1823
1824        while (!tstack.empty())
1825        {
1826                VspKdNode *node = tstack.top();
1827                tstack.pop();
1828
1829                if (node->IsLeaf())
1830                {
1831                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
1832
1833                        float c = leaf->GetAvgRayContribution();
1834                        int num = (int)(c*ratioPerLeaf + 0.5);
1835                        //                      cout<<num<<" ";
1836
1837                        for (int i=0; i < num; i++)
1838                        {
1839                                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1840                                /*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1841                                Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1842                                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1843                                rays.push_back(SimpleRay(origin, direction));*/
1844                        }
1845
1846                }
1847                else
1848                {
1849                        VspKdInterior *in =
1850                                dynamic_cast<VspKdInterior *>(node);
1851               
1852                        tstack.push(in->GetFront());
1853                        tstack.push(in->GetBack());
1854                }
1855        }
1856
1857        return (int)rays.size();
1858}
1859
1860
1861float VspKdTree::GetAvgPvsSize()
1862{
1863        stack<VspKdNode *> tstack;
1864        tstack.push(mRoot);
1865
1866        int sumPvs = 0;
1867        int leaves = 0;
1868
1869        while (!tstack.empty())
1870        {
1871                VspKdNode *node = tstack.top();
1872                tstack.pop();
1873
1874                if (node->IsLeaf())
1875                {
1876                        VspKdLeaf *leaf = (VspKdLeaf *)node;
1877                        // update pvs size
1878                        leaf->UpdatePvsSize();
1879                        sumPvs += leaf->GetPvsSize();
1880                        leaves++;
1881                }
1882                else
1883                {
1884                        VspKdInterior *in = (VspKdInterior *)node;
1885                       
1886                        tstack.push(in->GetFront());
1887                        tstack.push(in->GetBack());
1888                }
1889        }
1890
1891        return sumPvs / (float)leaves;
1892}
1893
1894VspKdNode *VspKdTree::GetRoot() const
1895{
1896        return mRoot;
1897}
1898
1899AxisAlignedBox3 VspKdTree::GetBBox(VspKdNode *node) const
1900{
1901        if (node->mParent == NULL)
1902                return mBox;
1903
1904        if (!node->IsLeaf())
1905        {
1906                if (node->Type() == VspKdNode::EIntermediate)
1907            return (dynamic_cast<VspKdInterior *>(node))->mBox;
1908                else
1909                        return (dynamic_cast<VspKdIntermediate *>(node))->mBox;
1910        }
1911
1912        VspKdInterior *parent = dynamic_cast<VspKdInterior *>(node->mParent);
1913
1914        if (parent->mAxis >= 3)
1915                return parent->mBox;
1916
1917        AxisAlignedBox3 box(parent->mBox);
1918        if (parent->GetFront() == node)
1919                box.SetMin(parent->mAxis, parent->mPosition);
1920        else
1921                box.SetMax(parent->mAxis, parent->mPosition);
1922
1923        return box;
1924}
1925
1926int     VspKdTree::GetRootPvsSize() const
1927{
1928        return GetPvsSize(mRoot, mBox);
1929}
1930
1931const VspKdStatistics &VspKdTree::GetStatistics() const
1932{
1933        return mStat;
1934}
1935
1936void VspKdTree::AddRays(VssRayContainer &add)
1937{
1938        VssRayContainer remove;
1939        UpdateRays(remove, add);
1940}
1941
1942// return memory usage in MB
1943float VspKdTree::GetMemUsage() const
1944{
1945        return
1946                (sizeof(VspKdTree) +
1947                 mStat.Leaves() * sizeof(VspKdLeaf) +
1948                 mStat.Interior() * sizeof(VspKdInterior) +
1949                 mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
1950}
1951
1952float VspKdTree::GetRayMemUsage() const
1953{
1954        return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
1955}
1956
1957
1958int VspKdTree::ComputePvsSize(VspKdNode *node,
1959                                                          const RayInfoContainer &globalRays) const
1960{
1961        int pvsSize = 0;
1962
1963        const AxisAlignedBox3 box = GetBBox(node);
1964
1965        RayInfoContainer::const_iterator it, it_end = globalRays.end();
1966
1967        Intersectable::NewMail();
1968
1969        // warning: implicit conversion from VssRay to Ray
1970        for (it = globalRays.begin(); it != globalRays.end(); ++ it)
1971        {
1972                VssRay *ray = (*it).mRay;
1973
1974                // ray intersects node bounding box
1975                if (box.GetMinMaxT(*ray, NULL, NULL))
1976                {
1977                        if (ray->mTerminationObject && !ray->mTerminationObject->Mailed())
1978                        {
1979                                ray->mTerminationObject->Mail();
1980                                ++ pvsSize;
1981                        }
1982                }
1983        }
1984        return pvsSize;
1985}
1986
1987
1988void VspKdTree::CollectLeaves(vector<VspKdLeaf *> &leaves) const
1989{
1990        stack<VspKdNode *> nodeStack;
1991        nodeStack.push(mRoot);
1992
1993        while (!nodeStack.empty())
1994        {
1995                VspKdNode *node = nodeStack.top();
1996
1997                nodeStack.pop();
1998
1999                if (node->IsLeaf())
2000                {
2001                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2002                        leaves.push_back(leaf);
2003                }
2004                else
2005                {
2006                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2007
2008                        nodeStack.push(interior->GetBack());
2009                        nodeStack.push(interior->GetFront());
2010                }
2011        }
2012}
2013
2014
2015int VspKdTree::FindNeighbors(VspKdLeaf *n,
2016                                                         vector<VspKdLeaf *> &neighbors,
2017                                                         bool onlyUnmailed)
2018{
2019        stack<VspKdNode *> nodeStack;
2020
2021        nodeStack.push(mRoot);
2022
2023        AxisAlignedBox3 box = GetBBox(n);
2024
2025        while (!nodeStack.empty())
2026        {
2027                VspKdNode *node = nodeStack.top();
2028                nodeStack.pop();
2029
2030                if (node->IsLeaf())
2031                {
2032                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2033
2034                        if (leaf != n && (!onlyUnmailed || !leaf->Mailed()))
2035                                neighbors.push_back(leaf);
2036                }
2037                else
2038                {
2039                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2040
2041                        if (interior->mPosition > box.Max(interior->mAxis))
2042                                nodeStack.push(interior->mBack);
2043                        else
2044                        {
2045                                if (interior->mPosition < box.Min(interior->mAxis))
2046                                        nodeStack.push(interior->mFront);
2047                                else
2048                                {
2049                                        // random decision
2050                                        nodeStack.push(interior->mBack);
2051                                        nodeStack.push(interior->mFront);
2052                                }
2053                        }
2054                }
2055        }
2056
2057        return (int)neighbors.size();
2058}
2059
2060
2061void VspKdTree::SetViewCellsManager(ViewCellsManager *vcm)
2062{
2063        mViewCellsManager = vcm;
2064}
2065
2066
2067int VspKdTree::CastLineSegment(const Vector3 &origin,
2068                                                           const Vector3 &termination,
2069                                                           vector<ViewCell *> &viewcells)
2070{
2071        int hits = 0;
2072
2073        float mint = 0.0f, maxt = 1.0f;
2074        const Vector3 dir = termination - origin;
2075
2076        stack<LineTraversalData> tStack;
2077
2078        Intersectable::NewMail();
2079
2080        Vector3 entp = origin;
2081        Vector3 extp = termination;
2082
2083        VspKdNode *node = mRoot;
2084        VspKdNode *farChild;
2085
2086        float position;
2087        int axis;
2088
2089        while (1)
2090        {
2091                if (!node->IsLeaf())
2092                {
2093                        VspKdInterior *in = dynamic_cast<VspKdInterior *>(node);
2094                        position = in->mPosition;
2095                        axis = in->mAxis;
2096
2097                        if (entp[axis] <= position)
2098                        {
2099                                if (extp[axis] <= position)
2100                                {
2101                                        node = in->mBack;
2102                                        // cases N1,N2,N3,P5,Z2,Z3
2103                                        continue;
2104                                } else
2105                                {
2106                                        // case N4
2107                                        node = in->mBack;
2108                                        farChild = in->mFront;
2109                                }
2110                        }
2111                        else
2112                        {
2113                                if (position <= extp[axis])
2114                                {
2115                                        node = in->mFront;
2116                                        // cases P1,P2,P3,N5,Z1
2117                                        continue;
2118                                }
2119                                else
2120                                {
2121                                        node = in->mFront;
2122                                        farChild = in->mBack;
2123                                        // case P4
2124                                }
2125                        }
2126
2127                        // $$ modification 3.5.2004 - hints from Kamil Ghais
2128                        // case N4 or P4
2129                        float tdist = (position - origin[axis]) / dir[axis];
2130                        tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
2131                        extp = origin + dir * tdist;
2132                        maxt = tdist;
2133                }
2134                else
2135                {
2136                        // compute intersection with all objects in this leaf
2137                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2138                        ViewCell *vc = leaf->GetViewCell();
2139
2140                        if (!vc->Mailed())
2141                        {
2142                                vc->Mail();
2143                                viewcells.push_back(vc);
2144                                ++ hits;
2145                        }
2146
2147                        if (0) //matt TODO: REMOVE LATER
2148                                leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
2149
2150                        // get the next node from the stack
2151                        if (tStack.empty())
2152                                break;
2153
2154                        entp = extp;
2155                        mint = maxt;
2156                       
2157                        LineTraversalData &s  = tStack.top();
2158                        node = s.mNode;
2159                        extp = s.mExitPoint;
2160                        maxt = s.mMaxT;
2161                        tStack.pop();
2162                }
2163        }
2164
2165        return hits;
2166}
2167
2168
2169void VspKdTree::GenerateViewCell(VspKdLeaf *leaf)
2170{
2171        RayInfoContainer::const_iterator it, it_end = leaf->GetRays().end();
2172
2173        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(mViewCellsManager->GenerateViewCell());
2174        leaf->SetViewCell(vc);
2175
2176        vc->SetVolume(GetBBox(leaf).GetVolume());
2177        vc->SetArea(GetBBox(leaf).SurfaceArea());
2178
2179        vc->mLeaf = leaf;
2180
2181        for (it = leaf->GetRays().begin(); it != it_end; ++ it)
2182        {
2183                VssRay *ray = (*it).mRay;
2184
2185                if (ray->mTerminationObject)
2186                  vc->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf);
2187
2188                if (ray->mOriginObject)
2189                  vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf);
2190        }
2191}
2192
2193
2194bool VspKdTree::MergeViewCells(VspKdLeaf *l1, VspKdLeaf *l2)
2195{
2196        //-- change pointer to view cells of all leaves associated
2197        //-- with the previous view cells
2198        VspKdViewCell *fVc = l1->GetViewCell();
2199        VspKdViewCell *bVc = l2->GetViewCell();
2200
2201        VspKdViewCell *vc = dynamic_cast<VspKdViewCell *>(
2202                mViewCellsManager->MergeViewCells(fVc, bVc));
2203
2204        // if merge was unsuccessful
2205        if (!vc) return false;
2206// TODO
2207#if HAS_TO_BE_REDONE
2208        // set new size of view cell
2209        vc->SetVolume(fVc->GetVolume() + bVc->GetVolume());
2210        // new area
2211        vc->SetArea(fVc->GetArea() + bVc->GetArea());
2212
2213        vector<VspKdLeaf *> fLeaves = fVc->mLeaves;
2214        vector<VspKdLeaf *> bLeaves = bVc->mLeaves;
2215
2216        vector<VspKdLeaf *>::const_iterator it;
2217
2218        //-- change view cells of all the other leaves the view cell belongs to
2219        for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
2220        {
2221                (*it)->SetViewCell(vc);
2222                vc->mLeaves.push_back(*it);
2223        }
2224
2225        for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
2226        {
2227                (*it)->SetViewCell(vc);
2228                vc->mLeaves.push_back(*it);
2229        }
2230#endif
2231        vc->Mail();
2232
2233        // clean up old view cells
2234        DEL_PTR(fVc);
2235        DEL_PTR(bVc);
2236
2237        return true;
2238}
2239
2240void VspKdTree::CollectMergeCandidates()
2241{
2242        VspKdMergeCandidate::sOverallCost = 0;
2243        vector<VspKdLeaf *> leaves;
2244
2245        // collect the leaves, e.g., the "voxels" that will build the view cells
2246        CollectLeaves(leaves);
2247
2248        VspKdLeaf::NewMail();
2249
2250        vector<VspKdLeaf *>::const_iterator it, it_end = leaves.end();
2251
2252        // generate view cells
2253        for (it = leaves.begin(); it != it_end; ++ it)
2254                GenerateViewCell(*it);
2255
2256        // find merge candidates and push them into queue
2257        for (it = leaves.begin(); it != it_end; ++ it)
2258        {
2259                VspKdLeaf *leaf = *it;
2260                ViewCell *vc = leaf->GetViewCell();
2261                // no leaf is part of two merge candidates
2262                leaf->Mail();
2263                VspKdMergeCandidate::sOverallCost +=
2264                        vc->GetVolume() * vc->GetPvs().GetSize();
2265                vector<VspKdLeaf *> neighbors;
2266                FindNeighbors(leaf, neighbors, true);
2267
2268                vector<VspKdLeaf *>::const_iterator nit,
2269                        nit_end = neighbors.end();
2270
2271                for (nit = neighbors.begin(); nit != nit_end; ++ nit)
2272                {
2273                        mMergeQueue.push(VspKdMergeCandidate(leaf, *nit));
2274                }
2275        }
2276}
2277
2278
2279int VspKdTree::MergeViewCells(const VssRayContainer &rays)
2280{
2281        CollectMergeCandidates();
2282
2283        int merged = 0;
2284
2285        int vcSize = mStat.Leaves();
2286        // use priority queue to merge leaves
2287        while (!mMergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
2288                   (mMergeQueue.top().GetMergeCost() <
2289                    mMergeMaxCostRatio * VspKdMergeCandidate::sOverallCost))
2290        {
2291                //Debug << "mergecost: " << mergeQueue.top().GetMergeCost() /
2292                //VspKdMergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
2293                VspKdMergeCandidate mc = mMergeQueue.top();
2294                mMergeQueue.pop();
2295
2296                // both view cells equal!
2297                if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
2298                        continue;
2299
2300                if (mc.Valid())
2301                {
2302                        ViewCell::NewMail();
2303                        MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
2304
2305                        ++ merged;
2306                        -- vcSize;
2307                        // increase absolute merge cost
2308                        VspKdMergeCandidate::sOverallCost += mc.GetMergeCost();
2309                }
2310                // merge candidate not valid, because one of the leaves was already
2311                // merged with another one
2312                else
2313                {
2314                        // validate and reinsert into queue
2315                        mc.SetValid();
2316                        mMergeQueue.push(mc);
2317                        //Debug << "validate " << mc.GetMergeCost() << endl;
2318                }
2319        }
2320
2321        Debug << "merged " << merged << " of " << mStat.Leaves() << " leaves" << endl;
2322
2323        //TODO: should return sample contributions
2324        return merged;
2325}
2326
2327
2328void VspKdTree::RepairViewCellsLeafLists()
2329{
2330        // list not valid anymore => clear
2331        stack<VspKdNode *> nodeStack;
2332        nodeStack.push(mRoot);
2333
2334        ViewCell::NewMail();
2335
2336        while (!nodeStack.empty())
2337        {
2338                VspKdNode *node = nodeStack.top();
2339                nodeStack.pop();
2340
2341                if (node->IsLeaf())
2342                {
2343                        VspKdLeaf *leaf = dynamic_cast<VspKdLeaf *>(node);
2344
2345                        VspKdViewCell *viewCell = leaf->GetViewCell();
2346#if HAS_TO_BE_REDONE
2347                        if (!viewCell->Mailed())
2348                        {
2349                                viewCell->mLeaves.clear();
2350                                viewCell->Mail();
2351                        }
2352
2353                        viewCell->mLeaves.push_back(leaf);
2354#endif
2355                }
2356                else
2357                {
2358                        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2359
2360                        nodeStack.push(interior->GetFront());
2361                        nodeStack.push(interior->GetBack());
2362                }
2363        }
2364}
2365
2366
2367VspKdNode *VspKdTree::CollapseTree(VspKdNode *node, int &collapsed)
2368{
2369        if (node->IsLeaf())
2370                return node;
2371
2372        VspKdInterior *interior = dynamic_cast<VspKdInterior *>(node);
2373
2374        VspKdNode *front = CollapseTree(interior->GetFront(), collapsed);
2375        VspKdNode *back = CollapseTree(interior->GetBack(), collapsed);
2376
2377        if (front->IsLeaf() && back->IsLeaf())
2378        {
2379                VspKdLeaf *frontLeaf = dynamic_cast<VspKdLeaf *>(front);
2380                VspKdLeaf *backLeaf = dynamic_cast<VspKdLeaf *>(back);
2381
2382                //-- collapse tree
2383                if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
2384                {
2385                        VspKdViewCell *vc = frontLeaf->GetViewCell();
2386
2387                        VspKdLeaf *leaf = new VspKdLeaf(interior->GetParent(), 0);
2388                        leaf->SetViewCell(vc);
2389                        //leaf->SetTreeValid(frontLeaf->TreeValid());
2390
2391                        // replace a link from node's parent
2392                        if (leaf->mParent)
2393                        {
2394                                // replace a link from node's parent
2395                                VspKdInterior *parent =
2396                                        dynamic_cast<VspKdInterior *>(leaf->mParent);
2397                                parent->ReplaceChildLink(node, leaf);
2398                        }
2399                        ++ collapsed;
2400                        delete interior;
2401
2402                        return leaf;
2403                }
2404        }
2405
2406        return node;
2407}
2408
2409
2410int VspKdTree::CollapseTree()
2411{
2412        int collapsed = 0;
2413        CollapseTree(mRoot, collapsed);
2414        // revalidate leaves
2415        RepairViewCellsLeafLists();
2416
2417        return collapsed;
2418}
2419
2420
2421int VspKdTree::RefineViewCells(const VssRayContainer &rays)
2422{
2423        int shuffled = 0;
2424
2425        Debug << "refining " << (int)mMergeQueue.size() << " candidates " << endl;
2426        BspLeaf::NewMail();
2427
2428        // Use priority queue of remaining leaf pairs
2429        // These candidates either share the same view cells or
2430        // are border leaves which share a boundary.
2431        // We test if they can be shuffled, i.e.,
2432        // either one leaf is made part of one view cell or the other
2433        // leaf is made part of the other view cell. It is tested if the
2434        // remaining view cells are "better" than the old ones.
2435        while (!mMergeQueue.empty())
2436        {
2437                VspKdMergeCandidate mc = mMergeQueue.top();
2438                mMergeQueue.pop();
2439
2440                // both view cells equal or already shuffled
2441                if ((mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell()) ||
2442                        (mc.GetLeaf1()->Mailed()) || (mc.GetLeaf2()->Mailed()))
2443                        continue;
2444               
2445                // candidate for shuffling
2446                const bool wasShuffled = false;
2447                        //ShuffleLeaves(mc.GetLeaf1(), mc.GetLeaf2());
2448                        // matt: restore
2449                //-- stats
2450                if (wasShuffled)
2451                        ++ shuffled;
2452        }
2453
2454        return shuffled;
2455}
2456
2457
2458inline int AddedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2459{
2460        return pvs1.AddPvs(pvs2);
2461}
2462
2463
2464inline int SubtractedPvsSize(ObjectPvs pvs1, const ObjectPvs &pvs2)
2465{
2466        return pvs1.SubtractPvs(pvs2);
2467}
2468
2469
2470float EvalShuffleCost(VspKdLeaf *leaf, VspKdViewCell *vc1, VspKdViewCell *vc2)
2471{
2472        const int pvs1 = SubtractedPvsSize(vc1->GetPvs(), *leaf->mPvs);
2473        const int pvs2 = AddedPvsSize(vc2->GetPvs(), *leaf->mPvs);
2474
2475        const float vol1 = vc1->GetVolume() - leaf->mVolume;
2476        const float vol2 = vc2->GetVolume() + leaf->mVolume;
2477
2478        const float cost1 = pvs1 * vol1;
2479        const float cost2 = pvs2 * vol2;
2480
2481        return cost1 + cost2;
2482}
2483
2484
2485void VspKdTree::ShuffleLeaf(VspKdLeaf *leaf,
2486                                                        VspKdViewCell *vc1,
2487                                                        VspKdViewCell *vc2) const
2488{
2489        // compute new pvs and area
2490#if HAS_TO_BE_REDONE
2491        vc1->GetPvs().SubtractPvs(*leaf->mPvs);
2492        vc2->GetPvs().AddPvs(*leaf->mPvs);
2493       
2494        vc1->SetArea(vc1->GetVolume() - leaf->mVolume);
2495        vc2->SetArea(vc2->GetVolume() + leaf->mVolume);
2496
2497        /// add to second view cell
2498        vc2->mLeaves.push_back(leaf);
2499
2500        // erase leaf from old view cell
2501        vector<VspKdLeaf *>::iterator it = vc1->mLeaves.begin();
2502
2503        for (; *it != leaf; ++ it);
2504        vc1->mLeaves.erase(it);
2505
2506        leaf->SetViewCell(vc2);  // finally change view cell
2507#endif
2508}
2509
2510
2511bool VspKdTree::ShuffleLeaves(VspKdLeaf *leaf1, VspKdLeaf *leaf2) const
2512{
2513        VspKdViewCell *vc1 = leaf1->GetViewCell();
2514        VspKdViewCell *vc2 = leaf2->GetViewCell();
2515
2516        const float cost1 = vc1->GetPvs().GetSize() * vc1->GetArea();
2517        const float cost2 = vc2->GetPvs().GetSize() * vc2->GetArea();
2518
2519        const float oldCost = cost1 + cost2;
2520       
2521        float shuffledCost1 = Limits::Infinity;
2522        float shuffledCost2 = Limits::Infinity;
2523
2524        // the view cell should not be empty after the shuffle
2525//todo
2526        /*      if (vc1->mLeaves.size() > 1)
2527                shuffledCost1 = EvalShuffleCost(leaf1, vc1, vc2);
2528        if (vc2->mLeaves.size() > 1)
2529                shuffledCost2 = EvalShuffleCost(leaf2, vc2, vc1);
2530
2531*/      // shuffling unsuccessful
2532        if ((oldCost <= shuffledCost1) && (oldCost <= shuffledCost2))
2533                return false;
2534       
2535        if (shuffledCost1 < shuffledCost2)
2536        {
2537                ShuffleLeaf(leaf1, vc1, vc2);
2538                leaf1->Mail();
2539        }
2540        else
2541        {
2542                ShuffleLeaf(leaf2, vc2, vc1);
2543                leaf2->Mail();
2544        }
2545
2546        return true;
2547}
2548
2549
2550
2551
2552/*********************************************************************/
2553/*                VspKdMergeCandidate implementation                      */
2554/*********************************************************************/
2555
2556
2557VspKdMergeCandidate::VspKdMergeCandidate(VspKdLeaf *l1, VspKdLeaf *l2):
2558mMergeCost(0),
2559mLeaf1(l1),
2560mLeaf2(l2),
2561mLeaf1Id(l1->GetViewCell()->mMailbox),
2562mLeaf2Id(l2->GetViewCell()->mMailbox)
2563{
2564        EvalMergeCost();
2565}
2566
2567float VspKdMergeCandidate::GetLeaf1Cost() const
2568{
2569        VspKdViewCell *vc = mLeaf1->GetViewCell();
2570        return vc->GetPvs().GetSize() * vc->GetVolume();
2571}
2572
2573float VspKdMergeCandidate::GetLeaf2Cost() const
2574{
2575        VspKdViewCell *vc = mLeaf2->GetViewCell();
2576        return vc->GetPvs().GetSize() * vc->GetVolume();
2577}
2578
2579void VspKdMergeCandidate::EvalMergeCost()
2580{
2581        //-- compute pvs difference
2582        VspKdViewCell *vc1 = mLeaf1->GetViewCell();
2583        VspKdViewCell *vc2 = mLeaf2->GetViewCell();
2584
2585        const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
2586        const int vcPvs = diff1 + vc1->GetPvs().GetSize();
2587
2588        //-- compute ratio of old cost
2589        //-- (added size of left and right view cell times pvs size)
2590        //-- to new rendering cost (size of merged view cell times pvs size)
2591        const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();
2592       
2593        const float newCost =
2594                (float)vcPvs * (vc1->GetVolume() + vc2->GetVolume());
2595
2596        mMergeCost = newCost - oldCost;
2597
2598//      if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
2599//              mMergeCost += 1.0;
2600}
2601
2602void VspKdMergeCandidate::SetLeaf1(VspKdLeaf *l)
2603{
2604        mLeaf1 = l;
2605}
2606
2607void VspKdMergeCandidate::SetLeaf2(VspKdLeaf *l)
2608{
2609        mLeaf2 = l;
2610}
2611
2612
2613VspKdLeaf *VspKdMergeCandidate::GetLeaf1()
2614{
2615        return mLeaf1;
2616}
2617
2618
2619VspKdLeaf *VspKdMergeCandidate::GetLeaf2()
2620{
2621        return mLeaf2;
2622}
2623
2624
2625bool VspKdMergeCandidate::Valid() const
2626{
2627        return
2628                (mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
2629                (mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
2630}
2631
2632
2633float VspKdMergeCandidate::GetMergeCost() const
2634{
2635        return mMergeCost;
2636}
2637
2638
2639void VspKdMergeCandidate::SetValid()
2640{
2641        mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
2642        mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;
2643
2644        EvalMergeCost();
2645}
Note: See TracBrowser for help on using the repository browser.