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

Revision 1359, 47.2 KB checked in by mattausch, 18 years ago (diff)

worked on global object sorting

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