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

Revision 425, 39.3 KB checked in by mattausch, 19 years ago (diff)
Line 
1
2// ================================================================
3// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4// ****************************************************************
5/**
6   The KD tree based LSDS
7*/
8// Initial coding by
9/**
10   @author Jiri Bittner
11*/
12
13// Standard headers
14#include <stack>
15#include <queue>
16#include <algorithm>
17#include <fstream>
18#include <string>
19
20#include "VspKdTree.h"
21
22#include "Environment.h"
23#include "VssRay.h"
24#include "Intersectable.h"
25#include "Ray.h"
26#include "RayInfo.h"
27
28// Static variables
29int VspKdTreeLeaf::mailID = 0;
30
31/// Adds object to the pvs of the front and back node
32inline void AddObject2Pvs(Intersectable *object,
33                                                  const int side,
34                                                  int &pvsBack,
35                                                  int &pvsFront)
36{
37        if (!object)
38                return;
39               
40        if (side <= 0)
41        {
42                if (!object->Mailed() && !object->Mailed(2))
43                {
44                        ++ pvsBack;
45
46                        if (object->Mailed(1))
47                                object->Mail(2);
48                        else
49                                object->Mail();
50                }
51        }
52       
53        if (side >= 0)
54        {
55                if (!object->Mailed(1) && !object->Mailed(2))
56                {
57                        ++ pvsFront;
58                        if (object->Mailed())
59                                object->Mail(2);
60                        else
61                                object->Mail(1);
62                }
63        }
64}
65
66/**************************************************************/
67/*                class VspKdTreeNode implementation          */
68/**************************************************************/
69
70// Inline constructor
71VspKdTreeNode::VspKdTreeNode(VspKdTreeInterior *p):
72mParent(p), mAxis(-1), mDepth(p ? p->mDepth + 1 : 0)
73{}
74
75VspKdTreeNode::~VspKdTreeNode()
76{};
77
78inline VspKdTreeInterior *VspKdTreeNode::GetParent() const
79{
80        return mParent;
81}
82
83inline void VspKdTreeNode::SetParent(VspKdTreeInterior *p)
84{
85        mParent = p;
86}
87
88bool VspKdTreeNode::IsLeaf() const
89{
90        return mAxis == -1;
91}
92 
93int VspKdTreeNode::GetAccessTime()
94{
95        return 0x7FFFFFF;
96}
97
98/**************************************************************/
99/*                VspKdTreeInterior implementation            */
100/**************************************************************/
101
102VspKdTreeInterior::VspKdTreeInterior(VspKdTreeInterior *p):
103VspKdTreeNode(p), mBack(NULL), mFront(NULL), mAccesses(0), mLastAccessTime(-1)
104{
105}
106
107int VspKdTreeInterior::GetAccessTime()
108{
109        return mLastAccessTime;
110}
111
112void VspKdTreeInterior::SetupChildLinks(VspKdTreeNode *b, VspKdTreeNode *f)
113{
114        mBack = b;
115        mFront = f;
116        b->SetParent(this);
117        f->SetParent(this);
118}
119
120void VspKdTreeInterior::ReplaceChildLink(VspKdTreeNode *oldChild,
121                                                                                 VspKdTreeNode *newChild)
122{
123        if (mBack == oldChild)
124                mBack = newChild;
125        else
126                mFront = newChild;
127}
128
129int VspKdTreeInterior::Type() const 
130{
131        return EInterior;
132}
133 
134VspKdTreeInterior::~VspKdTreeInterior()
135{
136        DEL_PTR(mBack);
137        DEL_PTR(mFront);
138}
139 
140void VspKdTreeInterior::Print(ostream &s) const
141{
142        switch (mAxis)
143        {
144        case 0:
145                s << "x "; break;
146        case 1:
147                s << "y "; break;
148        case 2:
149                s << "z "; break;
150        }
151
152        s << mPosition << " ";
153
154        mBack->Print(s);
155        mFront->Print(s);
156}
157       
158int VspKdTreeInterior::ComputeRayIntersection(const RayInfo &rayData, float &t)
159{
160        return rayData.ComputeRayIntersection(mAxis, mPosition, t);
161}
162
163VspKdTreeNode *VspKdTreeInterior::GetBack() const
164{
165        return mBack;
166}
167
168VspKdTreeNode *VspKdTreeInterior::GetFront() const
169{
170        return mFront;
171}
172
173
174/**************************************************************/
175/*              class VspKdTreeLeaf implementation            */
176/**************************************************************/
177VspKdTreeLeaf::VspKdTreeLeaf(VspKdTreeInterior *p,      const int nRays):
178VspKdTreeNode(p), mRays(), mPvsSize(0), mValidPvs(false)
179{
180        mRays.reserve(nRays);
181}
182
183VspKdTreeLeaf::~VspKdTreeLeaf()
184{}
185
186int VspKdTreeLeaf::Type() const 
187{
188        return ELeaf;
189}
190
191void VspKdTreeLeaf::Print(ostream &s) const
192{
193        s << endl << "L: r = " << (int)mRays.size() << endl;
194};
195
196void VspKdTreeLeaf::AddRay(const RayInfo &data)
197{
198        mValidPvs = false;
199        mRays.push_back(data);
200        data.mRay->Ref();
201}
202
203int VspKdTreeLeaf::GetPvsSize() const
204{
205        return mPvsSize;
206}
207
208void VspKdTreeLeaf::SetPvsSize(const int s)
209{
210        mPvsSize = s;
211}
212
213void VspKdTreeLeaf::Mail()
214{
215        mMailbox = mailID;
216}
217
218void VspKdTreeLeaf::NewMail()
219{
220        ++ mailID;
221}
222
223bool VspKdTreeLeaf::Mailed() const
224{
225        return mMailbox == mailID;
226}
227
228bool VspKdTreeLeaf::Mailed(const int mail)
229{
230        return mMailbox >= mailID + mail;
231}
232
233float VspKdTreeLeaf::GetAvgRayContribution() const
234{
235        return GetPvsSize() / ((float)mRays.size() + Limits::Small);
236}
237
238float VspKdTreeLeaf::GetSqrRayContribution() const
239{
240        return sqr(GetAvgRayContribution());
241}
242
243/*********************************************************/
244/*            class VspKdTree implementation             */
245/*********************************************************/
246
247
248// Constructor
249VspKdTree::VspKdTree():
250mOnlyDrivingAxis(true)
251{         
252        environment->GetIntValue("VspKdTree.Termination.maxDepth", mTermMaxDepth);
253        environment->GetIntValue("VspKdTree.Termination.minPvs", mTermMinPvs);
254        environment->GetIntValue("VspKdTree.Termination.minRays", mTermMinRays);
255        environment->GetFloatValue("VspKdTree.Termination.maxRayContribution", mTermMaxRayContribution);
256        environment->GetFloatValue("VspKdTree.Termination.maxCostRatio", mTermMaxCostRatio);
257        environment->GetFloatValue("VspKdTree.Termination.minSize", mTermMinSize);
258
259        mTermMinSize = sqr(mTermMinSize);
260
261        environment->GetFloatValue("VspKdTree.epsilon", mEpsilon);
262        environment->GetFloatValue("VspKdTree.ct_div_ci", mCtDivCi);
263       
264        environment->GetFloatValue("VspKdTree.maxTotalMemory", mMaxTotalMemory);
265        environment->GetFloatValue("VspKdTree.maxStaticMemory", mMaxStaticMemory);
266   
267        environment->GetIntValue("VspKdTree.accessTimeThreshold", mAccessTimeThreshold);
268        environment->GetIntValue("VspKdTree.minCollapseDepth", mMinCollapseDepth);
269
270        // split type
271        char sname[128];
272        environment->GetStringValue("VspKdTree.splitType", sname);
273        string name(sname);
274               
275        Debug << "======= vsp kd tree options ========" << endl;
276        Debug << "max depth: "<< mTermMaxDepth << endl;
277        Debug << "min pvs: "<< mTermMinPvs << endl;
278        Debug << "min rays: "<< mTermMinRays << endl;
279        Debug << "max ray contribution: "<< mTermMaxRayContribution << endl;
280        Debug << "max cost ratio: "<< mTermMaxCostRatio << endl;
281        Debug << "min size: "<<mTermMinSize << endl;
282
283        if (name.compare("regular") == 0)
284        {
285                Debug << "using regular split" << endl;
286                splitType = ESplitRegular;
287        }
288        else
289        {
290                if (name.compare("heuristic") == 0)
291                {
292                        Debug << "using heuristic split" << endl;
293                        splitType = ESplitHeuristic;
294                }
295                else
296                {
297                        cerr << "Invalid VspKdTree split type " << name << endl;
298                        exit(1);
299                }
300        }
301
302        mRoot = NULL;
303
304        mSplitCandidates = new vector<SortableEntry>;
305}
306
307
308VspKdTree::~VspKdTree()
309{
310        DEL_PTR(mRoot);
311}
312
313void VspKdStatistics::Print(ostream &app) const
314{
315  app << "===== VspKdTree statistics ===============\n";
316
317  app << "#N_RAYS ( Number of rays )\n"
318      << rays <<endl;
319 
320  app << "#N_INITPVS ( Initial PVS size )\n"
321      << initialPvsSize <<endl;
322 
323  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
324
325  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
326
327  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
328  for (int i=0; i<7; i++)
329    app << splits[i] <<" ";
330  app <<endl;
331
332  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
333    rayRefs << "\n";
334
335  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
336    rayRefs/(double)rays << "\n";
337
338  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
339    rayRefs/(double)Leaves() << "\n";
340
341  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
342    maxRayRefs << "\n";
343
344
345  //  app << setprecision(4);
346
347  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
348    maxDepthNodes*100/(double)Leaves()<<endl;
349
350  app << "#N_PMINPVSLEAVES  ( Percentage of leaves with minCost )\n"<<
351    minPvsNodes*100/(double)Leaves()<<endl;
352
353  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\n"<<
354    minRaysNodes*100/(double)Leaves()<<endl;
355       
356        app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
357    minSizeNodes*100/(double)Leaves()<<endl;
358
359        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
360    maxRayContribNodes*100/(double)Leaves()<<endl;
361
362  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
363    addedRayRefs<<endl;
364
365  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
366    removedRayRefs<<endl;
367
368  //  app << setprecision(4);
369
370  app << "#N_CTIME  ( Construction time [s] )\n"
371      << Time() << " \n";
372
373  app << "===== END OF VspKdTree statistics ==========\n";
374
375}
376
377
378void
379VspKdTreeLeaf::UpdatePvsSize()
380{
381        if (!mValidPvs)
382        {
383                Intersectable::NewMail();
384                int pvsSize = 0;
385                for(RayInfoContainer::iterator ri = mRays.begin();
386                        ri != mRays.end(); ++ ri)
387                {
388                        if ((*ri).mRay->IsActive())
389                        {
390                                Intersectable *object;
391#if BIDIRECTIONAL_RAY
392                                object = (*ri).mRay->mOriginObject;
393                       
394                                if (object && !object->Mailed())
395                                {
396                                        ++ pvsSize;
397                                        object->Mail();
398                                }
399#endif
400                                object = (*ri).mRay->mTerminationObject;
401                                if (object && !object->Mailed())
402                                {
403                                        ++ pvsSize;
404                                        object->Mail();
405                                }
406                        }
407                }
408
409                mPvsSize = pvsSize;
410                mValidPvs = true;
411        }
412}
413
414
415void
416VspKdTree::Construct(VssRayContainer &rays,
417                                         AxisAlignedBox3 *forcedBoundingBox)
418{
419        mStat.Start();
420 
421        mMaxMemory = mMaxStaticMemory;
422
423        DEL_PTR(mRoot);
424
425        mRoot = new VspKdTreeLeaf(NULL, (int)rays.size());
426
427        // first construct a leaf that will get subdivided
428        VspKdTreeLeaf *leaf = dynamic_cast<VspKdTreeLeaf *>(mRoot);
429       
430        mStat.nodes = 1;
431       
432        mBox.Initialize();
433       
434        for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
435        {
436                leaf->AddRay(RayInfo(*ri));
437
438                mBox.Include((*ri)->GetOrigin());
439                mBox.Include((*ri)->GetTermination());
440        }
441       
442        if (forcedBoundingBox)
443                mBox = *forcedBoundingBox;
444
445        Debug << "bbox = " << mBox << endl;
446       
447        mStat.rays = (int)leaf->mRays.size();
448        leaf->UpdatePvsSize();
449       
450        mStat.initialPvsSize = leaf->GetPvsSize();
451       
452        // Subdivide();
453        mRoot = Subdivide(TraversalData(leaf, mBox, 0));
454
455        if (mSplitCandidates)
456        {
457                // force release of this vector
458                delete mSplitCandidates;
459                mSplitCandidates = new vector<SortableEntry>;
460        }
461 
462        mStat.Stop();
463
464        mStat.Print(cout);
465        Debug << "#Total memory=" << GetMemUsage() << endl;
466}
467
468
469
470VspKdTreeNode *VspKdTree::Subdivide(const TraversalData &tdata)
471{
472        VspKdTreeNode *result = NULL;
473
474        priority_queue<TraversalData> tStack;
475        //stack<TraversalData> tStack;
476 
477        tStack.push(tdata);
478
479        AxisAlignedBox3 backBox;
480        AxisAlignedBox3 frontBox;
481 
482        int lastMem = 0;
483
484        while (!tStack.empty())
485        {
486                float mem = GetMemUsage();
487               
488                if (lastMem / 10 != ((int)mem) / 10)
489                {
490                        Debug << mem << " MB" << endl;
491                }
492                lastMem = (int)mem;
493               
494                if (1 && mem > mMaxMemory)
495                {
496                        Debug << "memory limit reached: " << mem << endl;
497                        // count statistics on unprocessed leafs
498                        while (!tStack.empty())
499                        {
500                                EvaluateLeafStats(tStack.top());
501                                tStack.pop();
502                        }
503                        break;
504                }
505   
506                TraversalData data = tStack.top();
507                tStack.pop();   
508               
509                VspKdTreeNode *node = SubdivideNode((VspKdTreeLeaf *) data.mNode,
510                                                                                        data.mBox, backBox,     frontBox);
511                if (result == NULL)
512                        result = node;
513   
514                if (!node->IsLeaf())
515                {
516                        VspKdTreeInterior *interior = dynamic_cast<VspKdTreeInterior *>(node);
517
518                        // push the children on the stack
519                        tStack.push(TraversalData(interior->GetBack(), backBox, data.mDepth + 1));
520                        tStack.push(TraversalData(interior->GetFront(), frontBox, data.mDepth + 1));
521                }
522                else
523                {
524                        EvaluateLeafStats(data);
525                }
526        }
527
528        return result;
529}
530
531
532// returns selected plane for subdivision
533int VspKdTree::SelectPlane(VspKdTreeLeaf *leaf,
534                                                   const AxisAlignedBox3 &box,
535                                                   float &position,
536                                                   int &raysBack,
537                                                   int &raysFront,
538                                                   int &pvsBack,
539                                                   int &pvsFront)
540{
541        int minDirDepth = 6;
542        int axis;
543        float costRatio;
544       
545        if (splitType == ESplitRegular)
546        {
547                costRatio = BestCostRatioRegular(leaf,
548                                                                                 axis,
549                                                                                 position,
550                                                                                 raysBack,
551                                                                                 raysFront,
552                                                                                 pvsBack,
553                                                                                 pvsFront);
554        }
555        else if (splitType == ESplitHeuristic)
556        {
557                        costRatio = BestCostRatioHeuristic(leaf,
558                                                                                           axis,       
559                                                                                           position,
560                                                                                           raysBack,
561                                                                                           raysFront,
562                                                                                           pvsBack,
563                                                                                           pvsFront);
564        }
565        else
566        {
567                cerr << "VspKdTree: Unknown split heuristics\n";
568                exit(1);
569        }
570
571        if (0 && costRatio > mTermMaxCostRatio)
572        {
573                Debug << "Too big cost ratio " << costRatio << endl;
574                return -1;
575        }
576
577#if 1   
578        Debug <<
579                "pvs=" << leaf->mPvsSize <<
580                " rays=" << (int)leaf->mRays.size() <<
581                " rc=" << leaf->GetAvgRayContribution() <<
582                " axis=" << axis << endl;
583#endif
584       
585        return axis;
586}
587
588
589                                                       
590float VspKdTree::EvalCostRatio(VspKdTreeLeaf *leaf,
591                                                           const int axis,
592                                                           const float position,
593                                                           int &raysBack,
594                                                           int &raysFront,
595                                                           int &pvsBack,
596                                                           int &pvsFront)
597{
598        raysBack = 0;
599        raysFront = 0;
600        pvsFront = 0;
601        pvsBack = 0;
602
603        float newCost;
604
605        Intersectable::NewMail(3);
606       
607        // eval pvs size
608        int pvsSize = leaf->GetPvsSize();
609
610        Intersectable::NewMail(3);
611
612        // this is the main ray classification loop!
613    for(RayInfoContainer::iterator ri = leaf->mRays.begin();
614                ri != leaf->mRays.end(); ++ ri)
615        {
616                if (!(*ri).mRay->IsActive())
617                        continue;
618                       
619                // determine the side of this ray with respect to the plane
620                int side = (*ri).ComputeRayIntersection(axis, position, (*ri).mRay->mT);
621                //                              (*ri).mRay->mSide = side;
622                       
623                if (side <= 0)
624                        ++ raysBack;
625                               
626                if (side >= 0)
627                        ++ raysFront;
628                               
629                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
630        }       
631
632        float ratio = 0;
633       
634        if (1)
635        {
636                const float sum = pvsBack + pvsFront;
637                const float oldCost = (float)pvsSize * 2;
638
639                return sum / oldCost;
640        }
641        else {
642                AxisAlignedBox3 box = GetBBox(leaf);
643       
644                float minBox = box.Min(axis);
645                float maxBox = box.Max(axis);
646                float sizeBox = maxBox - minBox;
647               
648                // float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
649                float sum = pvsBack * (position - minBox) + pvsFront * (maxBox - position);
650
651                newCost = mCtDivCi + sum  / sizeBox;
652       
653                //Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
654                //  float oldCost = leaf->mRays.size();
655                float oldCost = (float)pvsSize;
656               
657                float ratio = newCost / oldCost;
658       
659        }
660       
661        return ratio;
662}
663
664float VspKdTree::BestCostRatioRegular(VspKdTreeLeaf *leaf,
665                                                                          int &axis,
666                                                                          float &position,
667                                                                          int &raysBack,
668                                                                          int &raysFront,
669                                                                          int &pvsBack,
670                                                                          int &pvsFront)
671{
672        int nRaysBack[3], nRaysFront[3];
673        int nPvsBack[3], nPvsFront[3];
674
675        float nPosition[3];
676        float nCostRatio[3];
677        int bestAxis = -1;
678       
679        AxisAlignedBox3 sBox = GetBBox(leaf);
680
681        // int sAxis = box.Size().DrivingAxis();
682        int sAxis = sBox.Size().DrivingAxis();
683
684        for (axis = 0; axis < 3; ++ axis)
685        {
686                if (!mOnlyDrivingAxis || axis == sAxis)
687                {
688                        nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis]) * 0.5f;
689                                               
690                        nCostRatio[axis] = EvalCostRatio(leaf,
691                                                                                         axis,
692                                                                                         nPosition[axis],
693                                                                                         nRaysBack[axis],
694                                                                                         nRaysFront[axis],
695                                                                                         nPvsBack[axis],
696                                                                                         nPvsFront[axis]);
697                       
698                        if (bestAxis == -1)
699                                bestAxis = axis;
700                        else
701                                if (nCostRatio[axis] < nCostRatio[bestAxis])
702                                        bestAxis = axis;
703                }
704        }
705
706        axis = bestAxis;
707        position = nPosition[bestAxis];
708
709        raysBack = nRaysBack[bestAxis];
710        raysFront = nRaysFront[bestAxis];
711
712        pvsBack = nPvsBack[bestAxis];
713        pvsFront = nPvsFront[bestAxis];
714       
715        return nCostRatio[bestAxis];
716}
717
718float VspKdTree::BestCostRatioHeuristic(VspKdTreeLeaf *leaf,
719                                                                                int &axis,
720                                                                                float &position,
721                                                                                int &raysBack,
722                                                                                int &raysFront,
723                                                                                int &pvsBack,
724                                                                                int &pvsFront)
725{
726        AxisAlignedBox3 box = GetBBox(leaf);
727        //      AxisAlignedBox3 dirBox = GetDirBBox(node);
728       
729        axis = box.Size().DrivingAxis();
730               
731        SortSplitCandidates(leaf, axis);
732 
733        // go through the lists, count the number of objects left and right
734        // and evaluate the following cost funcion:
735        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
736       
737        int rl = 0, rr = (int)leaf->mRays.size();
738        int pl = 0, pr = leaf->GetPvsSize();
739
740        float minBox = box.Min(axis);
741        float maxBox = box.Max(axis);
742        float sizeBox = maxBox - minBox;
743       
744        float minBand = minBox + 0.1f*(maxBox - minBox);
745        float maxBand = minBox + 0.9f*(maxBox - minBox);
746       
747        float sum = rr*sizeBox;
748        float minSum = 1e20f;
749       
750        Intersectable::NewMail();
751
752        // set all object as belonging to the fron pvs
753        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
754                ri != leaf->mRays.end(); ++ ri)
755        {
756                if ((*ri).mRay->IsActive())
757                {
758                        Intersectable *object = (*ri).mRay->mTerminationObject;
759                       
760                        if (object)
761                                if (!object->Mailed())
762                                {
763                                        object->Mail();
764                                        object->mCounter = 1;
765                                }
766                                else
767                                        ++ object->mCounter;
768                }
769        }
770       
771        Intersectable::NewMail();
772       
773        for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
774                ci < mSplitCandidates->end(); ++ ci)
775        {
776                VssRay *ray;
777                       
778                switch ((*ci).type)
779                {
780                        case SortableEntry::ERayMin:
781                                {
782                                        ++ rl;
783                                        ray = (VssRay *) (*ci).data;
784                                        Intersectable *object = ray->mTerminationObject;
785                                        if (object && !object->Mailed())
786                                        {
787                                                object->Mail();
788                                                ++ pl;
789                                        }
790                                        break;
791                                }
792                        case SortableEntry::ERayMax:
793                                {
794                                        -- rr;
795                                       
796                                        ray = (VssRay *) (*ci).data;
797                                        Intersectable *object = ray->mTerminationObject;
798                                       
799                                        if (object)
800                                        {
801                                                if (-- object->mCounter == 0)
802                                                -- pr;
803                                        }
804                                               
805                                        break;
806                                }
807                }
808                       
809                if ((*ci).value > minBand && (*ci).value < maxBand)
810                {
811                        sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);
812               
813                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
814                        // cout<<"cost= "<<sum<<endl;
815                       
816                        if (sum < minSum)
817                        {
818                                minSum = sum;
819                                position = (*ci).value;
820                       
821                                raysBack = rl;
822                                raysFront = rr;
823                       
824                                pvsBack = pl;
825                                pvsFront = pr;
826                               
827                        }
828                }
829        }
830
831        float oldCost = (float)leaf->GetPvsSize();
832        float newCost = mCtDivCi + minSum / sizeBox;
833        float ratio = newCost / oldCost;
834 
835        //Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
836        //     <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
837       
838        return ratio;
839}
840
841void VspKdTree::SortSplitCandidates(VspKdTreeLeaf *node,
842                                                                        const int axis)
843{
844        mSplitCandidates->clear();
845 
846        int requestedSize = 2 * (int)(node->mRays.size());
847        // creates a sorted split candidates array
848        if (mSplitCandidates->capacity() > 500000 &&
849                requestedSize < (int)(mSplitCandidates->capacity()/10) )
850        {
851        delete mSplitCandidates;
852                mSplitCandidates = new vector<SortableEntry>;
853        }
854 
855        mSplitCandidates->reserve(requestedSize);
856
857        // insert all queries
858        for(RayInfoContainer::const_iterator ri = node->mRays.begin();
859                ri < node->mRays.end(); ++ ri)
860        {
861                bool positive = (*ri).mRay->HasPosDir(axis);
862                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
863                                                                                                  (*ri).ExtrapOrigin(axis), (void *)&*ri));
864               
865                mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
866                                                                                                  (*ri).ExtrapTermination(axis), (void *)&*ri));
867        }
868       
869        stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
870}
871
872
873void VspKdTree::EvaluateLeafStats(const TraversalData &data)
874{
875        // the node became a leaf -> evaluate stats for leafs
876        VspKdTreeLeaf *leaf = dynamic_cast<VspKdTreeLeaf *>(data.mNode);
877
878        if (data.mDepth >= mTermMaxDepth)
879                ++ mStat.maxDepthNodes;
880 
881        //  if ( (int)(leaf->mRays.size()) < termMinCost)
882        //    stat.minCostNodes++;
883        if (leaf->GetPvsSize() < mTermMinPvs)
884                ++ mStat.minPvsNodes;
885
886        if (leaf->GetPvsSize() < mTermMinRays)
887                ++ mStat.minRaysNodes;
888
889        if (0 && leaf->GetAvgRayContribution() > mTermMaxRayContribution)
890                ++ mStat.maxRayContribNodes;
891       
892        if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
893                ++ mStat.minSizeNodes;
894
895//      if ((int)(leaf->mRays.size()) > stat.maxRayRefs)
896//              mStat.maxRayRefs = (int)leaf->mRays.size();
897}
898
899
900inline bool VspKdTree::TerminationCriteriaMet(const VspKdTreeLeaf *leaf,
901                                                                                          const AxisAlignedBox3 &box) const
902{
903        return ((leaf->GetPvsSize() < mTermMinPvs) ||
904                    //(leaf->mRays.size() < mTermMinRays) ||
905                        // (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
906                        (leaf->mDepth >= mTermMaxDepth) ||
907                        (SqrMagnitude(box.Size()) <= mTermMinSize));
908}
909
910
911VspKdTreeNode *VspKdTree::SubdivideNode(VspKdTreeLeaf *leaf,
912                                                                                const AxisAlignedBox3 &box,
913                                                                                AxisAlignedBox3 &backBBox,
914                                                                                AxisAlignedBox3 &frontBBox)
915{
916        if (TerminationCriteriaMet(leaf, box))
917        {
918                if (1)
919                {
920                        if (leaf->mDepth >= mTermMaxDepth)
921                        {
922                                Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
923                                Debug << "Bbox: " << GetBBox(leaf) << endl;
924                        }
925               
926                        Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
927                }
928                return leaf;
929        }
930       
931        float position;
932        // first count ray sides
933        int raysBack;
934        int raysFront;
935        int pvsBack;
936        int pvsFront;
937       
938        // select subdivision axis
939        const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
940
941        Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" << pvsFront << endl;
942
943        if (axis == -1)
944        {
945                return leaf;
946        }
947 
948        mStat.nodes += 2;
949        //++ mStat.splits[axis];
950
951        // add the new nodes to the tree
952        VspKdTreeInterior *node = new VspKdTreeInterior(leaf->mParent);
953
954        node->mAxis = axis;
955        node->mPosition = position;
956        node->mBox = box;
957         
958        backBBox = box;
959        frontBBox = box;
960
961        VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
962        back->SetPvsSize(pvsBack);
963        VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
964        front->SetPvsSize(pvsFront);
965
966        // replace a link from node's parent
967        if (leaf->mParent)
968                leaf->mParent->ReplaceChildLink(leaf, node);
969        // and setup child links
970        node->SetupChildLinks(back, front);
971       
972        if (axis <= VspKdTreeNode::SPLIT_Z)
973        {
974                backBBox.SetMax(axis, position);
975                frontBBox.SetMin(axis, position);
976               
977                for(RayInfoContainer::iterator ri = leaf->mRays.begin();
978                                ri != leaf->mRays.end(); ++ ri)
979                {
980                        if ((*ri).mRay->IsActive())
981                        {
982                                // first unref ray from the former leaf
983                                (*ri).mRay->Unref();
984
985                                // determine the side of this ray with respect to the plane
986                                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
987
988                                if (side == 0)
989                                {
990                                        if ((*ri).mRay->HasPosDir(axis))
991                                        {
992                                                back->AddRay(RayInfo((*ri).mRay,
993                                                                                         (*ri).mMinT,
994                                                                                         (*ri).mRay->mT));
995                                                front->AddRay(RayInfo((*ri).mRay,
996                                                                                          (*ri).mRay->mT,
997                                                                                          (*ri).mMaxT));
998                                        }
999                                        else
1000                                        {
1001                                                back->AddRay(RayInfo((*ri).mRay,
1002                                                                                         (*ri).mRay->mT,
1003                                                                                         (*ri).mMaxT));
1004                                                front->AddRay(RayInfo((*ri).mRay,
1005                                                                                          (*ri).mMinT,
1006                                                                                          (*ri).mRay->mT));
1007                                        }
1008                                }
1009                                else
1010                                        if (side == 1)
1011                                                front->AddRay(*ri);
1012                                        else
1013                                                back->AddRay(*ri);
1014                        } else
1015                                (*ri).mRay->Unref();
1016                }
1017        }
1018        else
1019        {
1020                // rays front/back
1021   
1022        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1023                        ri != leaf->mRays.end(); ++ ri)
1024                {
1025                        if ((*ri).mRay->IsActive())
1026                        {
1027                                // first unref ray from the former leaf
1028                                (*ri).mRay->Unref();
1029
1030                                int side;
1031                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1032                                        side = 1;
1033                                else
1034                                        side = -1;
1035                               
1036                                if (side == 1)
1037                                        front->AddRay(*ri);
1038                                else
1039                                        back->AddRay(*ri);             
1040                        }
1041                        else
1042                                (*ri).mRay->Unref();
1043                }
1044        }
1045       
1046        // update stats
1047        mStat.rayRefs -= (int)leaf->mRays.size();
1048        mStat.rayRefs += raysBack + raysFront;
1049
1050        DEL_PTR(leaf);
1051
1052        return node;
1053}
1054
1055int VspKdTree::ReleaseMemory(const int time)
1056{
1057        stack<VspKdTreeNode *> tstack;
1058
1059        // find a node in the tree which subtree will be collapsed
1060        int maxAccessTime = time - mAccessTimeThreshold;
1061        int released = 0;
1062
1063        tstack.push(mRoot);
1064
1065        while (!tstack.empty())
1066        {
1067                VspKdTreeNode *node = tstack.top();
1068                tstack.pop();
1069   
1070                if (!node->IsLeaf())
1071                {
1072                        VspKdTreeInterior *in = dynamic_cast<VspKdTreeInterior *>(node);
1073                        // cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1074                        if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
1075                        {
1076                                released = CollapseSubtree(node, time);
1077                                break;
1078                        }
1079                        if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
1080                        {
1081                                tstack.push(in->GetFront());
1082                                tstack.push(in->GetBack());
1083                        }
1084                        else
1085                        {
1086                                tstack.push(in->GetBack());
1087                                tstack.push(in->GetFront());
1088                        }
1089                }
1090        }
1091
1092        while (tstack.empty())
1093        {
1094                // could find node to collaps...
1095                // cout<<"Could not find a node to release "<<endl;
1096                break;
1097        }
1098 
1099        return released;
1100}
1101
1102
1103VspKdTreeNode *VspKdTree::SubdivideLeaf(VspKdTreeLeaf *leaf,
1104                                                                                const float sizeThreshold)
1105{
1106        VspKdTreeNode *node = leaf;
1107
1108        AxisAlignedBox3 leafBBox = GetBBox(leaf);
1109
1110        static int pass = 0;
1111        ++ pass;
1112       
1113        // check if we should perform a dynamic subdivision of the leaf
1114        if (// leaf->mRays.size() > (unsigned)termMinCost &&
1115                (leaf->GetPvsSize() >= mTermMinPvs) &&
1116                (SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
1117        {
1118        // memory check and release
1119                if (GetMemUsage() > mMaxTotalMemory)
1120                        ReleaseMemory(pass);
1121               
1122                AxisAlignedBox3 backBBox, frontBBox;
1123
1124                // subdivide the node
1125                node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
1126        }
1127        return node;
1128}
1129
1130
1131
1132void VspKdTree::UpdateRays(VssRayContainer &remove,
1133                                                   VssRayContainer &add)
1134{
1135        VspKdTreeLeaf::NewMail();
1136
1137        // schedule rays for removal
1138        for(VssRayContainer::const_iterator ri = remove.begin();
1139                ri != remove.end(); ++ ri)
1140        {
1141                (*ri)->ScheduleForRemoval(); 
1142        }
1143
1144        int inactive = 0;
1145
1146        for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
1147        {
1148                if ((*ri)->ScheduledForRemoval())
1149                        //      RemoveRay(*ri, NULL, false);
1150                        // !!! BUG - with true it does not work correctly - aggreated delete
1151                        RemoveRay(*ri, NULL, true);
1152                else
1153                        ++ inactive;
1154        }
1155
1156        //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1157        for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
1158        {
1159                AddRay(*ri);
1160        }
1161}
1162
1163
1164void VspKdTree::RemoveRay(VssRay *ray,
1165                                                  vector<VspKdTreeLeaf *> *affectedLeaves,
1166                                                  const bool removeAllScheduledRays)
1167{
1168        stack<RayTraversalData> tstack;
1169
1170        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1171 
1172        RayTraversalData data;
1173
1174        // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1175        while (!tstack.empty())
1176        {
1177                data = tstack.top();
1178                tstack.pop();
1179
1180                if (!data.mNode->IsLeaf())
1181                {
1182                        // split the set of rays in two groups intersecting the
1183                        // two subtrees
1184                        TraverseInternalNode(data, tstack);
1185        }
1186                else
1187                {
1188                        // remove the ray from the leaf
1189                        // find the ray in the leaf and swap it with the last ray...
1190                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)data.mNode;
1191     
1192                        if (!leaf->Mailed())
1193                        {
1194                                leaf->Mail();
1195                               
1196                                if (affectedLeaves)
1197                                        affectedLeaves->push_back(leaf);
1198       
1199                                if (removeAllScheduledRays)
1200                                {
1201                                        int tail = (int)leaf->mRays.size() - 1;
1202
1203                                        for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
1204                                        {
1205                                                if (leaf->mRays[i].mRay->ScheduledForRemoval())
1206                                                {
1207                                                        // find a ray to replace it with
1208                                                        while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
1209                                                        {
1210                                                                ++ mStat.removedRayRefs;
1211                                                                leaf->mRays[tail].mRay->Unref();
1212                                                                leaf->mRays.pop_back();
1213                                                               
1214                                                                -- tail;
1215                                                        }
1216                                                       
1217                                                        if (tail < i)
1218                                                                break;
1219             
1220                                                        ++ mStat.removedRayRefs;
1221             
1222                                                        leaf->mRays[i].mRay->Unref();
1223                                                        leaf->mRays[i] = leaf->mRays[tail];
1224                                                        leaf->mRays.pop_back();
1225                                                        -- tail;
1226                                                }
1227                                        }
1228                                }
1229                        }
1230     
1231                        if (!removeAllScheduledRays)
1232                                for (int i=0; i < (int)leaf->mRays.size(); i++)
1233                                {
1234                                        if (leaf->mRays[i].mRay == ray)
1235                                        {
1236                                                ++ mStat.removedRayRefs;
1237                                                ray->Unref();
1238                                                leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
1239                                                leaf->mRays.pop_back();
1240                                            // check this ray again
1241                                                break;
1242                                        }
1243                                }
1244                }
1245        }
1246
1247        if (ray->RefCount() != 0)
1248        {
1249                cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1250                exit(1);
1251        }
1252
1253}
1254
1255void VspKdTree::AddRay(VssRay *ray)
1256{
1257        stack<RayTraversalData> tstack;
1258 
1259        tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1260 
1261        RayTraversalData data;
1262 
1263        while (!tstack.empty())
1264        {
1265                data = tstack.top();
1266                tstack.pop();
1267
1268                if (!data.mNode->IsLeaf())
1269                {
1270                        TraverseInternalNode(data, tstack);
1271                }
1272                else
1273                {
1274                        // remove the ray from the leaf
1275                        // find the ray in the leaf and swap it with the last ray
1276                        VspKdTreeLeaf *leaf = dynamic_cast<VspKdTreeLeaf *>(data.mNode);
1277
1278                        leaf->AddRay(data.mRayData);
1279                        ++ mStat.addedRayRefs;
1280                }
1281        }
1282}
1283
1284void VspKdTree::TraverseInternalNode(RayTraversalData &data,
1285                                                                         stack<RayTraversalData> &tstack)
1286{
1287        VspKdTreeInterior *in =  (VspKdTreeInterior *) data.mNode;
1288 
1289        if (in->mAxis <= VspKdTreeNode::SPLIT_Z)
1290        {
1291            // determine the side of this ray with respect to the plane
1292                int side = in->ComputeRayIntersection(data.mRayData, data.mRayData.mRay->mT);
1293   
1294      if (side == 0)
1295          {
1296                if (data.mRayData.mRay->HasPosDir(in->mAxis))
1297                {
1298                        tstack.push(RayTraversalData(in->GetBack(),
1299                                                RayInfo(data.mRayData.mRay, data.mRayData.mMinT, data.mRayData.mRay->mT)));
1300                               
1301                        tstack.push(RayTraversalData(in->GetFront(),
1302                                                RayInfo(data.mRayData.mRay, data.mRayData.mRay->mT, data.mRayData.mMaxT)));
1303       
1304                }
1305                else
1306                {
1307                        tstack.push(RayTraversalData(in->GetBack(),
1308                                                                                 RayInfo(data.mRayData.mRay,
1309                                                                                                 data.mRayData.mRay->mT,
1310                                                                                                 data.mRayData.mMaxT)));
1311                        tstack.push(RayTraversalData(in->GetFront(),
1312                                                                                 RayInfo(data.mRayData.mRay,
1313                                                                                                 data.mRayData.mMinT,
1314                                                                                                 data.mRayData.mRay->mT)));
1315                }
1316          }
1317          else
1318          if (side == 1)
1319                          tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1320                  else
1321                          tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1322        }
1323        else
1324        {
1325                // directional split
1326                if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
1327                        tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1328                else
1329                        tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1330        }
1331}
1332
1333
1334int VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
1335{
1336        // first count all rays in the subtree
1337        // use mail 1 for this purpose
1338        stack<VspKdTreeNode *> tstack;
1339       
1340        int rayCount = 0;
1341        int totalRayCount = 0;
1342        int collapsedNodes = 0;
1343
1344#if DEBUG_COLLAPSE
1345        cout << "Collapsing subtree" << endl;
1346        cout << "acessTime=" << sroot->GetAccessTime() << endl;
1347        cout << "depth=" << (int)sroot->depth << endl;
1348#endif
1349
1350        // tstat.collapsedSubtrees++;
1351        // tstat.collapseDepths += (int)sroot->depth;
1352        // tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1353        tstack.push(sroot);
1354        VssRay::NewMail();
1355       
1356        while (!tstack.empty())
1357        {
1358                ++ collapsedNodes;
1359               
1360                VspKdTreeNode *node = tstack.top();
1361                tstack.pop();
1362                if (node->IsLeaf())
1363                {
1364                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
1365                       
1366                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1367                                        ri != leaf->mRays.end(); ++ ri)
1368                        {
1369                                ++ totalRayCount;
1370                               
1371                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed())
1372                                {
1373                                        (*ri).mRay->Mail();
1374                                        ++ rayCount;
1375                                }
1376                        }
1377                }
1378                else
1379                {
1380                        tstack.push(((VspKdTreeInterior *)node)->GetFront());
1381                        tstack.push(((VspKdTreeInterior *)node)->GetBack());
1382                }
1383        }
1384 
1385        VssRay::NewMail();
1386
1387        // create a new node that will hold the rays
1388        VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf(sroot->mParent, rayCount);
1389       
1390        if (newLeaf->mParent)
1391                newLeaf->mParent->ReplaceChildLink(sroot, newLeaf);
1392
1393        tstack.push(sroot);
1394
1395        while (!tstack.empty())
1396        {
1397                VspKdTreeNode *node = tstack.top();
1398                tstack.pop();
1399
1400                if (node->IsLeaf())
1401                {
1402                        VspKdTreeLeaf *leaf = dynamic_cast<VspKdTreeLeaf *>(node);
1403
1404                        for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1405                                        ri != leaf->mRays.end(); ++ ri)
1406                        {
1407                                // unref this ray from the old node             
1408                                if ((*ri).mRay->IsActive())
1409                                {
1410                                        (*ri).mRay->Unref();
1411                                        if (!(*ri).mRay->Mailed())
1412                                        {
1413                                                (*ri).mRay->Mail();
1414                                                newLeaf->AddRay(*ri);
1415                                        }
1416                                }
1417                                else
1418                                        (*ri).mRay->Unref();
1419                        }
1420                }
1421                else
1422                {
1423                        VspKdTreeInterior *interior =
1424                                dynamic_cast<VspKdTreeInterior *>(node);
1425                        tstack.push(interior->GetBack());
1426                        tstack.push(interior->GetFront());
1427                }
1428        }
1429 
1430        // delete the node and all its children
1431        DEL_PTR(sroot);
1432 
1433        // for(VspKdTreeNode::SRayContainer::iterator ri = newleaf->mRays.begin();
1434    //      ri != newleaf->mRays.end(); ++ ri)
1435        //     (*ri).ray->UnMail(2);
1436
1437#if DEBUG_COLLAPSE
1438        cout << "Total memory before=" << GetMemUsage() << endl;
1439#endif
1440
1441        mStat.nodes -= collapsedNodes - 1;
1442        mStat.rayRefs -= totalRayCount - rayCount;
1443 
1444#if DEBUG_COLLAPSE
1445        cout << "collapsed nodes" << collapsedNodes << endl;
1446        cout << "collapsed rays" << totalRayCount - rayCount << endl;
1447        cout << "Total memory after=" << GetMemUsage() << endl;
1448        cout << "================================" << endl;
1449#endif
1450
1451        //  tstat.collapsedNodes += collapsedNodes;
1452        //  tstat.collapsedRays += totalRayCount - rayCount;
1453    return totalRayCount - rayCount;
1454}
1455
1456
1457int VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
1458{
1459        stack<VspKdTreeNode *> tstack;
1460        tstack.push(mRoot);
1461
1462        Intersectable::NewMail();
1463        int pvsSize = 0;
1464       
1465        while (!tstack.empty())
1466        {
1467                VspKdTreeNode *node = tstack.top();
1468                tstack.pop();
1469   
1470          if (node->IsLeaf())
1471                {
1472                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1473                        for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1474                                 ri != leaf->mRays.end(); ++ ri)
1475                        {
1476                                if ((*ri).mRay->IsActive())
1477                                {
1478                                        Intersectable *object;
1479#if BIDIRECTIONAL_RAY
1480                                        object = (*ri).mRay->mOriginObject;
1481                                       
1482                                        if (object && !object->Mailed())
1483                                        {
1484                                                ++ pvsSize;
1485                                                object->Mail();
1486                                        }
1487#endif
1488                                        object = (*ri).mRay->mTerminationObject;
1489                                        if (object && !object->Mailed())
1490                                        {
1491                                                ++ pvsSize;
1492                                                object->Mail();
1493                                        }
1494                                }
1495                        }
1496                }
1497                else
1498                {
1499                        VspKdTreeInterior *in = dynamic_cast<VspKdTreeInterior *>(node);
1500       
1501                        if (in->mAxis < 3)
1502                        {
1503                                if (box.Max(in->mAxis) >= in->mPosition)
1504                                        tstack.push(in->GetFront());
1505                               
1506                                if (box.Min(in->mAxis) <= in->mPosition)
1507                                        tstack.push(in->GetBack());
1508                        }
1509                        else
1510                        {
1511                                // both nodes for directional splits
1512                                tstack.push(in->GetFront());
1513                                tstack.push(in->GetBack());
1514                        }
1515                }
1516        }
1517
1518        return pvsSize;
1519}
1520
1521void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1522                                                                                         float &maxRayContribution,
1523                                                                                         float &avgRayContribution)
1524{
1525        stack<VspKdTreeNode *> tstack;
1526        tstack.push(mRoot);
1527       
1528        minRayContribution = 1.0f;
1529        maxRayContribution = 0.0f;
1530
1531        float sumRayContribution = 0.0f;
1532        int leaves = 0;
1533
1534        while (!tstack.empty())
1535        {
1536                VspKdTreeNode *node = tstack.top();
1537                tstack.pop();
1538                if (node->IsLeaf())
1539                {
1540                        leaves++;
1541                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1542                        float c = leaf->GetAvgRayContribution();
1543
1544                        if (c > maxRayContribution)
1545                                maxRayContribution = c;
1546                        if (c < minRayContribution)
1547                                minRayContribution = c;
1548                        sumRayContribution += c;
1549                       
1550                }
1551                else {
1552                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
1553                        // both nodes for directional splits
1554                        tstack.push(in->GetFront());
1555                        tstack.push(in->GetBack());
1556                }
1557        }
1558       
1559        Debug << "sum=" << sumRayContribution << endl;
1560        Debug << "leaves=" << leaves << endl;
1561        avgRayContribution = sumRayContribution / (float)leaves;
1562}
1563
1564
1565int VspKdTree::GenerateRays(const float ratioPerLeaf,
1566                                                        SimpleRayContainer &rays)
1567{
1568        stack<VspKdTreeNode *> tstack;
1569        tstack.push(mRoot);
1570       
1571        while (!tstack.empty())
1572        {
1573                VspKdTreeNode *node = tstack.top();
1574                tstack.pop();
1575               
1576                if (node->IsLeaf())
1577                {
1578                        VspKdTreeLeaf *leaf = dynamic_cast<VspKdTreeLeaf *>(node);
1579
1580                        float c = leaf->GetAvgRayContribution();
1581                        int num = (int)(c*ratioPerLeaf + 0.5);
1582                        //                      cout<<num<<" ";
1583
1584                        for (int i=0; i < num; i++)
1585                        {
1586                                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1587                                /*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1588                                Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1589                                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1590                                rays.push_back(SimpleRay(origin, direction));*/
1591                        }
1592                       
1593                }
1594                else
1595                {
1596                        VspKdTreeInterior *in =
1597                                dynamic_cast<VspKdTreeInterior *>(node);
1598                        // both nodes for directional splits
1599                        tstack.push(in->GetFront());
1600                        tstack.push(in->GetBack());
1601                }
1602        }
1603
1604        return (int)rays.size();
1605}
1606
1607
1608float VspKdTree::GetAvgPvsSize()
1609{
1610        stack<VspKdTreeNode *> tstack;
1611        tstack.push(mRoot);
1612
1613        int sumPvs = 0;
1614        int leaves = 0;
1615       
1616        while (!tstack.empty())
1617        {
1618                VspKdTreeNode *node = tstack.top();
1619                tstack.pop();
1620               
1621                if (node->IsLeaf())
1622                {
1623                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1624                        // update pvs size
1625                        leaf->UpdatePvsSize();
1626                        sumPvs += leaf->GetPvsSize();
1627                        leaves++;
1628                }
1629                else
1630                {
1631                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
1632                        // both nodes for directional splits
1633                        tstack.push(in->GetFront());
1634                        tstack.push(in->GetBack());
1635                }
1636        }
1637
1638        return sumPvs / (float)leaves;
1639}
1640
1641VspKdTreeNode *VspKdTree::GetRoot() const
1642{
1643        return mRoot;
1644}
1645
1646AxisAlignedBox3 VspKdTree::GetBBox(VspKdTreeNode *node) const
1647{
1648        if (node->mParent == NULL)
1649                return mBox;
1650
1651        if (!node->IsLeaf())
1652                return (dynamic_cast<VspKdTreeInterior *>(node))->mBox;
1653
1654        if (node->mParent->mAxis >= 3)
1655                return node->mParent->mBox;
1656   
1657        AxisAlignedBox3 box(node->mParent->mBox);
1658        if (node->mParent->mFront == node)
1659                box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
1660        else
1661                box.SetMax(node->mParent->mAxis, node->mParent->mPosition);
1662
1663        return box;
1664}
1665
1666int     VspKdTree::GetRootPvsSize() const
1667{
1668        return GetPvsSize(mRoot, mBox);
1669}
1670
1671const VspKdStatistics &VspKdTree::GetStatistics() const
1672{
1673        return mStat;
1674}
1675
1676void VspKdTree::AddRays(VssRayContainer &add)
1677{
1678        VssRayContainer remove;
1679        UpdateRays(remove, add);
1680}
1681 
1682// return memory usage in MB
1683float VspKdTree::GetMemUsage() const
1684{
1685        return
1686                (sizeof(VspKdTree) +
1687                 mStat.Leaves() * sizeof(VspKdTreeLeaf) +
1688                 mStat.Interior() * sizeof(VspKdTreeInterior) +
1689                 mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
1690}
1691       
1692float VspKdTree::GetRayMemUsage() const
1693{
1694        return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
1695}
1696
1697int VspKdTree::ComputePvsSize(VspKdTreeNode *node, const RayInfoContainer &globalRays) const
1698{
1699        int pvsSize = 0;
1700
1701        const AxisAlignedBox3 box = GetBBox(node);
1702
1703        RayInfoContainer::const_iterator it, it_end = globalRays.end();
1704
1705        for (it = globalRays.begin(); it != globalRays.end(); ++ it)
1706                ;//pvsSize += box.GetMinMaxT((*it).mRay, NULL, NULL);
1707
1708        return pvsSize;
1709}
Note: See TracBrowser for help on using the repository browser.