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

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