source: trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp @ 431

Revision 431, 41.3 KB checked in by mattausch, 19 years ago (diff)

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