source: GTP/trunk/Lib/Vis/Preprocessing/src/VssTree.cpp @ 1004

Revision 1004, 47.1 KB checked in by mattausch, 18 years ago (diff)

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