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

Revision 2105, 47.3 KB checked in by bittner, 18 years ago (diff)

simple ray separated

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