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

Revision 433, 40.8 KB checked in by mattausch, 19 years ago (diff)

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