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

Revision 863, 59.5 KB checked in by mattausch, 18 years ago (diff)

working on preprocessor integration
added iv stuff

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