source: GTP/trunk/Lib/Vis/Preprocessing/src/RssTree.cpp @ 2575

Revision 2575, 74.4 KB checked in by bittner, 17 years ago (diff)

big merge: preparation for havran ray caster, check if everything works

RevLine 
[447]1// ================================================================
2// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
3// ****************************************************************
4/**
5   The KD tree based LSDS
6*/
7// Initial coding by
8/**
9   @author Jiri Bittner
10*/
11
12// Standard headers
13#include <stack>
14#include <queue>
15#include <algorithm>
16#include <fstream>
17#include <string>
18
19#include "RssTree.h"
20
21#include "Environment.h"
22#include "VssRay.h"
23#include "Intersectable.h"
24#include "Ray.h"
25#include "Containers.h"
[492]26#include "ViewCell.h"
27#include "Exporter.h"
28#include "Preprocessor.h"
29#include "SceneGraph.h"
[1883]30#include "SamplingStrategy.h"
[1891]31#include "ViewCellsManager.h"
[447]32
[860]33
[863]34namespace GtpVisibilityPreprocessor {
[860]35
36
[447]37#define DEBUG_SPLIT_COST 0
38#define DEBUG_SPLITS 0
39
[1199]40#define EVAL_VIEWCELLS 0
41 
[447]42// Static variables
43int
44RssTreeLeaf::mailID = 0;
45
46inline void
47AddObject2Pvs(Intersectable *object,
48                          const int side,
49                          int &pvsBack,
50                          int &pvsFront)
51{
52 
53  if (!object)
54        return;
55 
56  if (side <= 0) {
57        if (!object->Mailed() && !object->Mailed(2)) {
58          pvsBack++;
59          if (object->Mailed(1))
60                object->Mail(2);
61          else
62                object->Mail();
63        }
64  }
65 
66  if (side >= 0) {
67        if (!object->Mailed(1) && !object->Mailed(2)) {
68          pvsFront++;
69          if (object->Mailed())
70                object->Mail(2);
71          else
72                object->Mail(1);
73        }
74  }
75}
76
[492]77inline void
78AddViewcells2Pvs(const ViewCellContainer &viewcells,
79                                 const int side,
80                                 int &viewcellsBack,
81                                 int &viewcellsFront)
82{
83  ViewCellContainer::const_iterator it = viewcells.begin();
84 
85  for (; it != viewcells.end(); ++it) {
86        ViewCell *viewcell = *it;
87        if (side <= 0) {
88          if (!viewcell->Mailed() && !viewcell->Mailed(2)) {
89                viewcellsBack++;
90                if (viewcell->Mailed(1))
91                  viewcell->Mail(2);
92                else
93                  viewcell->Mail();
94          }
95        }
96       
97        if (side >= 0) {
98          if (!viewcell->Mailed(1) && !viewcell->Mailed(2)) {
99                viewcellsFront++;
100                if (viewcell->Mailed())
101                  viewcell->Mail(2);
102                else
103                  viewcell->Mail(1);
104          }
105        }
106  }
107}
108
109
[447]110// Constructor
111RssTree::RssTree()
112{
[1004]113  Environment::GetSingleton()->GetIntValue("RssTree.maxDepth", termMaxDepth);
114  Environment::GetSingleton()->GetIntValue("RssTree.minPvs", termMinPvs);
115  Environment::GetSingleton()->GetIntValue("RssTree.minRays", termMinRays);
116  Environment::GetSingleton()->GetFloatValue("RssTree.maxRayContribution", termMaxRayContribution);
117  Environment::GetSingleton()->GetFloatValue("RssTree.maxCostRatio", termMaxCostRatio);
[447]118
[1004]119  Environment::GetSingleton()->GetFloatValue("RssTree.minSize", termMinSize);
[447]120  termMinSize = sqr(termMinSize);
121       
[1004]122  Environment::GetSingleton()->GetFloatValue("RssTree.refDirBoxMaxSize", refDirBoxMaxSize);
[447]123  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
124 
[1004]125  Environment::GetSingleton()->GetFloatValue("RssTree.epsilon", epsilon);
126  Environment::GetSingleton()->GetFloatValue("RssTree.ct_div_ci", ct_div_ci);
[447]127       
[1004]128  Environment::GetSingleton()->GetFloatValue("RssTree.maxTotalMemory", maxTotalMemory);
129  Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
[447]130 
[1004]131  Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
[447]132
133
134 
[1004]135  Environment::GetSingleton()->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
[447]136  //= 1000;
[1004]137  Environment::GetSingleton()->GetIntValue("RssTree.minCollapseDepth", minCollapseDepth);
[447]138  //  int minCollapseDepth = 4;
139
140  //  pRefDirThresh = cos(0.5*M_PI - M_PI*refDirAngle/180.0);
141  //  cosRefDir = cos(M_PI*refDirAngle/180.0);
142  //  sinRefDir = sin(M_PI*refDirAngle/180.0);
143 
144 
145  // split type
146  char sname[128];
[1004]147  Environment::GetSingleton()->GetStringValue("RssTree.splitType", sname);
[447]148  string name(sname);
149       
150  if (name.compare("regular") == 0)
151    splitType = ESplitRegular;
152  else
153    if (name.compare("heuristic") == 0)
154      splitType = ESplitHeuristic;
155        else
156          if (name.compare("hybrid") == 0)
157                splitType = ESplitHybrid;
158          else {
159                cerr<<"Invalid RssTree split type "<<name<<endl;
160                exit(1);
161          }
[1199]162
163  Environment::GetSingleton()->GetIntValue("RssTree.hybridDepth", mHybridDepth);
164
[1004]165  Environment::GetSingleton()->GetBoolValue("RssTree.randomize", randomize);
166  Environment::GetSingleton()->GetBoolValue("RssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
[447]167
[1004]168  Environment::GetSingleton()->GetBoolValue("RssTree.interleaveDirSplits", mInterleaveDirSplits);
169  Environment::GetSingleton()->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);
[447]170
[1004]171  Environment::GetSingleton()->GetBoolValue("RssTree.importanceBasedCost", mImportanceBasedCost);
[492]172 
[1004]173  Environment::GetSingleton()->GetIntValue("RssTree.maxRays", mMaxRays);
[492]174
[1004]175  Environment::GetSingleton()->GetBoolValue("RssTree.perObjectTree", mPerObjectTree);
[447]176 
[516]177  //  mRoots;
178 
[447]179  splitCandidates = new vector<SortableEntry>;
180}
181
182
183RssTree::~RssTree()
184{
[516]185  for (int i=0; i < mRoots.size(); i++)
186        if (mRoots[i])
187          delete mRoots[i];
[447]188}
189
190
191
192
193void
194RssStatistics::Print(ostream &app) const
195{
[492]196  app << "###### RssTree statistics ######\n";
[447]197
198  app << "#N_RAYS ( Number of rays )\n"
199      << rays <<endl;
200
201  app << "#N_INITPVS ( Initial PVS size )\n"
202      << initialPvsSize <<endl;
203 
204  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
205
206  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
207
[667]208  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
[447]209  for (int i=0; i<7; i++)
210    app << splits[i] <<" ";
211  app <<endl;
212
213  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
214    rayRefs << "\n";
215
216  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
217    rayRefs/(double)rays << "\n";
218
219  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
220    rayRefs/(double)Leaves() << "\n";
221
222  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
223    maxRayRefs << "\n";
224
225
226  //  app << setprecision(4);
227
228  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
229    maxDepthNodes*100/(double)Leaves()<<endl;
230
[551]231  app << "#N_PMINPVSLEAVES  ( Percentage of leaves with minPvs )\n"<<
[447]232    minPvsNodes*100/(double)Leaves()<<endl;
233
[551]234  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minRays )\n"<<
[447]235    minRaysNodes*100/(double)Leaves()<<endl;
236       
237  app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
238    minSizeNodes*100/(double)Leaves()<<endl;
239
240  app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
241    maxRayContribNodes*100/(double)Leaves()<<endl;
242
243  app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
244    maxCostRatioNodes*100/(double)Leaves()<<endl;
245
246  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
247    addedRayRefs<<endl;
248
249  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
250    removedRayRefs<<endl;
251
252  //  app << setprecision(4);
253
254  app << "#N_CTIME  ( Construction time [s] )\n"
255      << Time() << " \n";
256
[492]257  app << "###### END OF RssTree statistics ######\n";
[447]258
259}
260
261
262void
[492]263RssTree::UpdatePvsSize(RssTreeLeaf *leaf)
[447]264{
[492]265  if (!leaf->mValidPvs) {
[447]266        Intersectable::NewMail();
267        int pvsSize = 0;
[492]268        for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
269                ri != leaf->rays.end();
[447]270                ri++)
271          if ((*ri).mRay->IsActive()) {
272                Intersectable *object;
[492]273                object = (*ri).GetObject();
[447]274                if (object && !object->Mailed()) {
275                  pvsSize++;
276                  object->Mail();
277                }
278          }
[492]279        leaf->SetPvsSize(pvsSize);
280       
281        ComputeImportance(leaf);
[447]282  }
283}
284
285bool
286RssTree::ClipRay(
287                                 RssTreeNode::RayInfo &rayInfo,
288                                 const AxisAlignedBox3 &box
289                                 )
290{
291  float tmin, tmax;
292  static Ray ray;
[492]293  cerr<<"Clip not reimplented yet...Quiting\n";
294  exit(1);
295  ray.Init(rayInfo.GetOrigin(), rayInfo.GetDir(), Ray::LINE_SEGMENT);
296
297  if (!box.ComputeMinMaxT(ray, &tmin, &tmax) || tmin>=tmax)
[447]298        return false;
299 
300  // now check if the ray origin lies inside the box
301  if ( tmax < rayInfo.mRay->GetSize() ) {
302        // this ray does not leave the box
303        rayInfo.mRay->SetupEndPoints(
304                                                                 ray.Extrap(tmax),
305                                                                 rayInfo.mRay->mTermination
306                                                                 );
307        return true;
308  }
309
310  return false;
311}
312
313
314void
315RssTree::Construct(
[516]316                                   ObjectContainer &objects,
317                                   VssRayContainer &rays
[447]318                                   // forced bounding box is only used when computing from-box
319                                   // visibility
[516]320                                   //                              AxisAlignedBox3 *forcedBoundingBox
[447]321                                   )
322{
[516]323  cout<<"Constructing rss tree"<<endl<<flush;
[447]324  stat.Start();
325 
326  maxMemory = maxStaticMemory;
327
[516]328  //  if (root)
329  //    delete root;
[492]330
[516]331  if (mPerObjectTree) {
332        // get max id from the rays
333        int i;
334        mRoots.resize(objects.size());
335        for (i = 0; i < objects.size(); i++) {
336          RssTreeLeaf *leaf = new RssTreeLeaf(NULL, 0);
337          //      leaf->bbox.Initialize();
338          leaf->dirBBox.Initialize();
339          leaf->dirBBox.SetMin(2, 0.0f);
[1824]340          leaf->dirBBox.SetMax(2, 0.0f);
[516]341          mRoots[i] = leaf;
342        }
[1112]343
[516]344        // init spatial bounding boxes
345        for (i = 0; i < objects.size(); i++) {
346          GetRoot(objects[i])->bbox = objects[i]->GetBox();
347        }
348        stat.nodes = i;
[569]349        stat.leaves = i;
[516]350  } else {
351        mRoots.resize(1);
[1563]352        RssTreeLeaf *leaf = new RssTreeLeaf(NULL, (int)rays.size());
[516]353        leaf->bbox.Initialize();
354        leaf->dirBBox.Initialize();
355        leaf->dirBBox.SetMin(2, 0.0f);
[1867]356        leaf->dirBBox.SetMax(2, 0.0f);
[516]357        mRoots[0] = leaf;
358        stat.nodes = 1;
[569]359        stat.leaves = 1;
[516]360  }
[447]361
362  for(VssRayContainer::const_iterator ri = rays.begin();
[516]363          ri != rays.end();
364          ri++) {
365       
366          RssTreeNode::RayInfo info(*ri);
[1900]367          if (!EpsilonEqualV3(info.GetOrigin(), info.GetTermination(), Limits::Small)) {
368
[516]369          // first construct a leaf that will get subdivide
[569]370          RssTreeLeaf *leaf = (RssTreeLeaf *) GetRoot(info.GetSourceObject());
[1900]371         
[447]372
[516]373          leaf->AddRay(info);
374         
375          // leaf bbox contains bbox of origins only
[1112]376
377          // include both origin and terminatin in the global bbox
378#if USE_ORIGIN
[516]379          leaf->bbox.Include((*ri)->GetOrigin());
380          bbox.Include((*ri)->GetOrigin());
[1112]381#else
[516]382          bbox.Include((*ri)->GetTermination());
[1112]383          leaf->bbox.Include((*ri)->GetTermination());
384#endif   
[516]385          Vector3 dVec = Vector3(
386                                                         (*ri)->GetDirParametrization(0),
387                                                         (*ri)->GetDirParametrization(1),
388                                                         0
389                                                         );
390
391          leaf->dirBBox.Include(dVec);
392          dirBBox.Include(dVec);
[1900]393          }
[447]394  }
[516]395
[492]396  // make the z axis (unused) a unit size
397  // important for volume computation
398 
[516]399  //  if ( forcedBoundingBox )
400  //    bbox = *forcedBoundingBox;
401 
[447]402  cout<<"Bbox = "<<bbox<<endl;
403  cout<<"Dirr Bbox = "<<dirBBox<<endl;
[516]404 
[2130]405  stat.rays = (int)rays.size();
[516]406  stat.initialPvsSize = 0;
407  for (int i=0; i < mRoots.size(); i++) {
408        RssTreeLeaf *leaf = (RssTreeLeaf *)mRoots[i];
[1112]409        UpdatePvsSize(leaf);
[516]410        stat.initialPvsSize += leaf->GetPvsSize();
411        mRoots[i] = Subdivide(TraversalData(leaf, GetBBox(leaf), 0));
412  }
413 
[447]414  if (splitCandidates) {
415    // force realease of this vector
416    delete splitCandidates;
417    splitCandidates = new vector<SortableEntry>;
418  }
419 
420  stat.Stop();
421  stat.Print(cout);
422  cout<<"#Total memory="<<GetMemUsage()<<endl;
423
[492]424  // this rotine also updates importances etc...
[447]425}
426
427int
428RssTree::UpdateSubdivision()
429{
430  priority_queue<TraversalData> tStack;
431  //  stack<TraversalData> tStack;
432 
[516]433  //  tStack.push(TraversalData(root, bbox, 0));
434  PushRoots(tStack);
435                       
[447]436  AxisAlignedBox3 backBox;
437  AxisAlignedBox3 frontBox;
438
439  maxMemory = maxTotalMemory;
440  int subdivided = 0;
441  int lastMem = 0;
442  while (!tStack.empty()) {
443               
444        float mem = GetMemUsage();
445               
446        if ( lastMem/10 != ((int)mem)/10) {
447          cout<<mem<<" MB"<<endl;
448        }
449        lastMem = (int)mem;
450               
451        if (  mem > maxMemory ) {
452      // count statistics on unprocessed leafs
453      while (!tStack.empty()) {
454                //                              EvaluateLeafStats(tStack.top());
455                tStack.pop();
456      }
457      break;
458    }
459   
460    TraversalData data = tStack.top();
461    tStack.pop();
462
463        if (data.node->IsLeaf()) {
464          RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
465                                                                                data.bbox,
466                                                                                backBox,
467                                                                                frontBox
468                                                                                );
469          if (!node->IsLeaf()) {
470                subdivided++;
[492]471               
[1112]472               
[447]473                RssTreeInterior *interior = (RssTreeInterior *) node;
474                // push the children on the stack
475                tStack.push(TraversalData(interior->back, backBox, data.depth+1));
476                tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
477          } else {
478                //      EvaluateLeafStats(data);
479          }
480        } else {
481          RssTreeInterior *interior = (RssTreeInterior *) data.node;
482          tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
483          tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
484        }
485  }
486  return subdivided;
487}
488
489
490RssTreeNode *
491RssTree::Subdivide(const TraversalData &tdata)
492{
493  RssTreeNode *result = NULL;
494
495  priority_queue<TraversalData> tStack;
496  //  stack<TraversalData> tStack;
497 
498  tStack.push(tdata);
499
500  AxisAlignedBox3 backBox;
501  AxisAlignedBox3 frontBox;
502
503 
504  int lastMem = 0;
505  while (!tStack.empty()) {
506
507        float mem = GetMemUsage();
508               
509        if ( lastMem/10 != ((int)mem)/10) {
510          cout<<mem<<" MB"<<endl;
511        }
512        lastMem = (int)mem;
513               
514        if (  mem > maxMemory ) {
515      // count statistics on unprocessed leafs
516      while (!tStack.empty()) {
517                EvaluateLeafStats(tStack.top());
518                tStack.pop();
519      }
520      break;
521    }
522   
523    TraversalData data = tStack.top();
524    tStack.pop();
525
526#if DEBUG_SPLITS
527        Debug<<"#Splitting node"<<endl;
528        data.node->Print(Debug);
529#endif
530        RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
531                                                                          data.bbox,
532                                                                          backBox,
533                                                                          frontBox
534                                                                          );
535       
536
537        if (result == NULL)
538      result = node;
539   
540    if (!node->IsLeaf()) {
541                       
542      RssTreeInterior *interior = (RssTreeInterior *) node;
543      // push the children on the stack
544      tStack.push(TraversalData(interior->back, backBox, data.depth+1));
545      tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
546
547#if DEBUG_SPLITS
548          Debug<<"#New nodes"<<endl;
549          interior->back->Print(Debug);
550          interior->front->Print(Debug);
551          Debug<<"#####################################"<<endl;
552#endif
553         
554    } else {
555      EvaluateLeafStats(data);
556    }
557  }
558
559  return result;
560}
561
562
563// returns selected plane for subdivision
564int
565RssTree::SelectPlane(
566                                         RssTreeLeaf *leaf,
567                                         const AxisAlignedBox3 &box,
[492]568                                         SplitInfo &info
[447]569                                         )
570{
571
[492]572  BestCostRatio(leaf,
573                                info);
574 
[447]575#if DEBUG_SPLIT_COST
[492]576  Debug<<"Split Info:"<<endl;
577  Debug<<"axis="<<info.axis<<" ratio="<<info.costRatio<<endl;
[1199]578  Debug<<"rays="<<info.rays<<
579        " rays back="<<info.raysBack<<
580        " rays front="<<info.raysFront<<endl;
581  Debug<<"c="<<info.contribution<<
582        " c back="<<info.contributionBack<<
583        " c front="<<info.contributionFront<<endl;
584
585  //  Debug<<"viewcells="<<info.viewCells<<
586  //    " viewcells back="<<info.viewCellsBack<<
587  //    " viewcells back="<<info.viewCellsFront<<endl;
[447]588#endif
[492]589 
590  if (info.costRatio > termMaxCostRatio) {
[447]591        //              cout<<"Too big cost ratio "<<costRatio<<endl;
592        stat.maxCostRatioNodes++;
593        return -1;
594  }
595       
[1199]596#if 0
[447]597  cout<<
598        "pvs="<<leaf->mPvsSize<<
599        " rays="<<leaf->rays.size()<<
600        " rc="<<leaf->GetAvgRayContribution()<<
[492]601        " axis="<<info.axis<<endl;
[447]602#endif
[492]603 
604  return info.axis;
[447]605}
606
607
[492]608void
[447]609RssTree::GetCostRatio(
610                                          RssTreeLeaf *leaf,
[492]611                                          SplitInfo &info
[447]612                                          )
613{
614
615  AxisAlignedBox3 box;
616  float minBox, maxBox;
617
[492]618  if (info.axis < 3) {
[447]619        box = GetBBox(leaf);
[492]620        minBox = box.Min(info.axis);
621        maxBox = box.Max(info.axis);
[447]622  }     else {
623        box = GetDirBBox(leaf);
[492]624        minBox = box.Min(info.axis-3);
625        maxBox = box.Max(info.axis-3);
[447]626  }
627       
628  float sizeBox = maxBox - minBox;
629
630  int pvsSize = leaf->GetPvsSize();
631
[492]632  if (!mImportanceBasedCost) {
633        const int costMethod = 0;
634       
635        switch (costMethod) {
636        case 0: {
637          float sum = info.pvsBack*(info.position - minBox) + info.pvsFront*(maxBox - info.position);
638          float newCost = ct_div_ci + sum/sizeBox;
[1500]639          float oldCost = (float)pvsSize;
[492]640          info.costRatio = newCost/oldCost;
641          break;
642        }
[1199]643        case 6: {
644          float sum = info.raysBack*(info.position - minBox) + info.raysFront*(maxBox - info.position);
645          float newCost = ct_div_ci + sum/sizeBox;
[1548]646          float oldCost = (float)info.rays;
[1199]647          info.costRatio = newCost/oldCost;
648          break;
649        }
650        case 7: {
651          float benefit =
652                (sqr(info.contributionBack) + sqr(info.contributionFront))/
653                (2.0f*sqr(info.contribution));
654          info.costRatio = 1.0f/(benefit + Limits::Small);
655          break;
656        }
657
658         
[492]659        case 1: {
660          float newContrib =
[1112]661                (info.contributionBack*(info.position - minBox) +
662                 +
663                 info.contributionFront*(maxBox - info.position))/sizeBox;
664          float oldContrib = info.contribution;
665          info.costRatio = newContrib/oldContrib;
666          //      cout<<info.contribution<<" "<<info.contributionBack<<" "<<info.contributionFront<<endl;
[492]667          break;
668        }
669        case 2: {
670          float sum =
671                info.viewCellsBack*(info.position - minBox) +
672                info.viewCellsFront*(maxBox - info.position);
673          float newCost = ct_div_ci + sum/sizeBox;
[1548]674          float oldCost = (float)info.viewCells;
[492]675          info.costRatio = newCost/oldCost;
676          break;
677        }
678        case 3: {
[1548]679          float newCost = (float)(info.raysBack*info.pvsBack  + info.raysFront*info.pvsFront);
[1189]680          float oldCost = (float)leaf->rays.size()*pvsSize;
[492]681          info.costRatio = newCost/oldCost;
[1877]682          break;
[492]683        }
684        }
[447]685  } else {
[492]686        const int costMethod = 1;
[447]687        // importance based cost
[492]688        switch (costMethod) {
689        case 0: {
690          break;
691        }
692        case 1: {
693          float newContrib =
694                sqr(info.pvsBack/(info.raysBack + Limits::Small)) +
695                sqr(info.pvsFront/(info.raysFront + Limits::Small));
696          float oldContrib = sqr(leaf->GetAvgRayContribution());
697          info.costRatio = oldContrib/newContrib;
698          break;
699        }
700        case 2: {
701          float newCost = (info.pvsBack  + info.pvsFront)*0.5f;
[1548]702          float oldCost = (float)pvsSize;
[492]703          info.costRatio = newCost/oldCost;
704          break;
705        }
706        case 3: {
[1548]707          float newCost = (float)abs(info.raysBack  - info.raysFront);
[1189]708          float oldCost = (float)leaf->rays.size();
[492]709          info.costRatio = newCost/oldCost;
710          break;
711        }
[447]712  }
[492]713  }
[447]714}
715                                                       
716
[492]717void
[447]718RssTree::EvalCostRatio(
719                                           RssTreeLeaf *leaf,
[492]720                                           SplitInfo &info
[447]721                                           )
722{
[1199]723  info.rays = 0;
[492]724  info.raysBack = 0;
725  info.raysFront = 0;
726  info.pvsFront = 0;
727  info.pvsBack = 0;
728  info.viewCells = 0;
729  info.viewCellsFront = 0;
730  info.viewCellsBack = 0;
[447]731
732
[1199]733
[492]734  float sumWeights = Limits::Small;
735  float sumWeightsBack = Limits::Small;
736  float sumWeightsFront = Limits::Small;
737 
738  float sumContribution = 0.0f;
739  float sumContributionBack = 0.0f;
740  float sumContributionFront = 0.0f;
[447]741
[492]742  Intersectable::NewMail();
[447]743
[1199]744#if EVAL_VIEWCELLS
[492]745  // count the numebr of viewcells first
746  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
747          ri != leaf->rays.end();
748          ri++)
749        if ((*ri).mRay->IsActive()) {
[1199]750
[492]751          ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
752          for (; it != (*ri).mRay->mViewCells.end(); ++it) {
753                if (!(*it)->Mailed()) {
754                  (*it)->Mail();
755                  info.viewCells++;
756                }
757          }
758        }
[1199]759#endif
[492]760 
761  Intersectable::NewMail(3);
762 
763  // this is the main ray classification loop!
764  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
765          ri != leaf->rays.end();
766          ri++)
767        if ((*ri).mRay->IsActive()) {
768          int side;
769         
770         
771          // determine the side of this ray with respect to the plane
772          side = (*ri).ComputeRaySide(info.axis, info.position);
773         
774          float weight, contribution;
775         
776          GetRayContribution(*ri, weight, contribution);
777          sumWeights += weight;
778          sumContribution += contribution;
[1199]779          info.rays++;
[492]780          if (side <= 0) {
781                info.raysBack++;
782                sumWeightsBack += weight;
783                sumContributionBack += contribution;
784          }
785         
786          if (side >= 0) {
787                info.raysFront++;
788                sumWeightsFront += weight;
789                sumContributionFront += contribution;
790          }
791         
792          AddObject2Pvs((*ri).GetObject(), side, info.pvsBack, info.pvsFront);
[1199]793
794#if EVAL_VIEWCELLS
[492]795          AddViewcells2Pvs((*ri).mRay->mViewCells,
796                                           side,
797                                           info.viewCellsBack,
798                                           info.viewCellsFront);
[1199]799#endif
[492]800        }
801
[1112]802  if (sumWeights!=0.0f)
803        info.contribution = sumContribution/sumWeights;
804  else
805        info.contribution = 0.0f;
806
807  if (sumWeightsBack!=0.0f)
808        info.contributionBack = sumContributionBack/sumWeightsBack;
809  else
810        info.contributionBack = 0.0f;
811
812  if (sumWeightsFront!=0.0f)
813        info.contributionFront = sumContributionFront/sumWeightsFront;
814  else
815        info.contributionFront = 0.0f;
[1199]816
817  //  info.costRatio = 0.1 + Random(0.5f);
818  //  return;
819
[492]820  GetCostRatio(
821                           leaf,
822                           info);
823 
[1112]824  //cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
[447]825  //  float oldCost = leaf->rays.size();
[492]826 
[1112]827  //  cout<<"ratio="<<ratio<<endl;
[447]828}
829
[492]830void
[447]831RssTree::BestCostRatio(
832                                           RssTreeLeaf *leaf,
[492]833                                           SplitInfo &info
[447]834                                           )
835{
[492]836  SplitInfo nInfo[5];
[447]837  int bestAxis = -1;
[492]838 
[447]839  AxisAlignedBox3 sBox = GetBBox(leaf);
840  AxisAlignedBox3 dBox = GetDirBBox(leaf);
841  // int sAxis = box.Size().DrivingAxis();
842  int sAxis = sBox.Size().DrivingAxis();
843  int dAxis = dBox.Size().DrivingAxis() + 3;
[492]844 
[447]845  float dirSplitBoxSize = 0.01f;
846  bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
847               
[1199]848  bool useCyclingAxis = true;
849
850  if (splitType == ESplitHeuristic ||
851          (splitType == ESplitHybrid && leaf->depth > mHybridDepth))
852        useCyclingAxis = false;
[447]853 
[1199]854  int cyclingAxis = 0;
855 
856  if (leaf->parent)
857        cyclingAxis = (leaf->parent->axis + 1)%5;
858 
859  if (0 && useCyclingAxis) {
860        int axis = cyclingAxis;
861        info.position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
862        info.axis = axis;
863        info.costRatio = 0.5f; // not true but good for speeding up the subdivision at the top levels!
[2130]864        info.raysBack = info.raysFront = (int)leaf->rays.size()/2;
[1199]865        return;
866  }
867
868 
[492]869  for (int axis = 0; axis < 5; axis++)
[447]870        if (
871                (axis < 3 && (leaf->depth < mDirSplitDepth ||  mInterleaveDirSplits)) ||
872                (axis >= 3 && (leaf->depth >= mDirSplitDepth))
873                ) {
[492]874          nInfo[axis].axis = axis;
[1199]875          if ((!useCyclingAxis && (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis))
876                  || (useCyclingAxis && (axis == cyclingAxis))
877                  ) {
[447]878               
879                if (splitType == ESplitRegular) {
880                  if (axis < 3)
[492]881                        nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
[447]882                  else
[492]883                        nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
884                  EvalCostRatio(leaf,
885                                                nInfo[axis]);
[447]886                } else
887                  if (splitType == ESplitHeuristic) {
[492]888                        EvalCostRatioHeuristic(
889                                                                   leaf,
890                                                                   nInfo[axis]
891                                                                   );
[447]892                  } else
893                        if (splitType == ESplitHybrid) {
[1199]894                          if (leaf->depth > mHybridDepth)
[492]895                                EvalCostRatioHeuristic(
896                                                                           leaf,
897                                                                           nInfo[axis]
898                                                                           );
[447]899                          else {
900                                if (axis < 3)
[492]901                                  nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
[447]902                                else
[492]903                                  nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
[447]904                               
[492]905                                EvalCostRatio(leaf,
906                                                          nInfo[axis]
907                                                          );
[447]908                          }
909                        } else {
910                          cerr<<"RssTree: Unknown split heuristics\n";
911                          exit(1);
912                        }
913               
914                if ( bestAxis == -1)
915                  bestAxis = axis;
916                else
[492]917                  if ( nInfo[axis].costRatio < nInfo[bestAxis].costRatio )
[447]918                        bestAxis = axis;
919          }
920        }
921 
[492]922  info = nInfo[bestAxis];
[447]923}
924
925       
[492]926void
[447]927RssTree::EvalCostRatioHeuristic(
928                                                                RssTreeLeaf *leaf,
[492]929                                                                SplitInfo &info
[447]930                                                                )
931{
932  AxisAlignedBox3 box;
933  float minBox, maxBox;
934       
[492]935  if (info.axis < 3) {
[447]936        box = GetBBox(leaf);
[492]937        minBox = box.Min(info.axis);
938        maxBox = box.Max(info.axis);
[447]939  } else {
940        box = GetDirBBox(leaf);
[492]941        minBox = box.Min(info.axis-3);
942        maxBox = box.Max(info.axis-3);
[447]943  }
944       
[1233]945  SortSubdivisionCandidates(leaf, info.axis);
[447]946 
947  // go through the lists, count the number of objects left and right
948  // and evaluate the following cost funcion:
949  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
[492]950 
951  SplitInfo currInfo;
952  currInfo.axis = info.axis;
953
954  currInfo.raysBack = 0;
[2130]955  currInfo.raysFront = (int)leaf->rays.size();
[492]956 
957  currInfo.pvsBack = 0;
958  currInfo.pvsFront = leaf->GetPvsSize();
959
[1199]960
961  float sumWeights = Limits::Small;
962  float sumWeightsBack = Limits::Small;
963  float sumWeightsFront = Limits::Small;
964 
965  float sumContribution = 0.0f;
966  float sumContributionBack = 0.0f;
967  float sumContributionFront = 0.0f;
968
969 
[447]970  float sizeBox = maxBox - minBox;
[492]971 
[485]972  float minBand = minBox + 0.1f*(maxBox - minBox);
973  float maxBand = minBox + 0.9f*(maxBox - minBox);
[447]974       
[492]975  // best cost ratio
976  info.costRatio = 1e20f;
977
978  currInfo.viewCells = 0;
[1199]979  currInfo.rays = 0;
[492]980
[447]981  Intersectable::NewMail();
982  // set all object as belonging to the fron pvs
983  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
984          ri != leaf->rays.end();
985          ri++)
986        if ((*ri).mRay->IsActive()) {
[1199]987
988          float weight, contribution;
989          GetRayContribution(*ri, weight, contribution);
990
991          sumWeights += weight;
992          sumContribution += contribution;
993          currInfo.rays++;
994
[492]995          Intersectable *object = (*ri).GetObject();
[447]996          if (object)
997                if (!object->Mailed()) {
998                  object->Mail();
[2053]999                 
[447]1000                  object->mCounter = 1;
1001                } else
[2053]1002                {
1003                       
[447]1004                  object->mCounter++;
[2053]1005                }
[1199]1006#if EVAL_VIEWCELLS       
[492]1007          // and do the same for all viewcells
1008          ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
1009         
1010          for (; it != (*ri).mRay->mViewCells.end(); ++it) {
1011                ViewCell *viewcell = *it;
1012                if (!viewcell->Mailed()) {
1013                  currInfo.viewCells++;
1014                  viewcell->Mail();
[2053]1015                 
[492]1016                  viewcell->mCounter = 1;
1017                } else
[2053]1018                {
1019                       
[492]1020                  viewcell->mCounter++;
[2053]1021                }
[492]1022          }
[1199]1023#endif
[447]1024        }
[492]1025
[1199]1026  if (sumWeights!=0.0f)
1027        currInfo.contribution = sumContribution/sumWeights;
1028  else
1029        currInfo.contribution = 0.0f;
1030
1031  sumWeightsFront = sumWeights;
1032  sumContributionFront = sumContribution;
1033
1034  sumWeightsBack = 0;
1035  sumContributionBack = 0;
1036 
[492]1037  currInfo.viewCellsBack = 0;
1038  currInfo.viewCellsFront = currInfo.viewCells;
[1199]1039 
[447]1040  Intersectable::NewMail();
[1199]1041 
[447]1042  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
1043          ci < splitCandidates->end();
1044          ci++) {
1045        switch ((*ci).type) {
1046        case SortableEntry::ERayMin: {
[492]1047          currInfo.raysFront--;
1048          currInfo.raysBack++;
[1199]1049
[492]1050          RssTreeNode::RayInfo *rayInfo = (RssTreeNode::RayInfo *) (*ci).data;
[1199]1051         
1052          float weight, contribution;
1053          GetRayContribution(*rayInfo, weight, contribution);
1054
1055          sumWeightsBack += weight;
1056          sumContributionBack += contribution;
1057
1058          sumWeightsFront -= weight;
1059          sumContributionFront -= contribution;
1060
[492]1061          Intersectable *object = rayInfo->GetObject();
[447]1062          if (object) {
1063                if (!object->Mailed()) {
1064                  object->Mail();
[492]1065                  currInfo.pvsBack++;
[447]1066                }
1067                if (--object->mCounter == 0)
[2053]1068                {
[492]1069                  currInfo.pvsFront--;
[2053]1070                }
[447]1071          }
[1199]1072
1073#if EVAL_VIEWCELLS       
[492]1074          ViewCellContainer::const_iterator it = rayInfo->mRay->mViewCells.begin();
1075          for (; it != rayInfo->mRay->mViewCells.end(); ++it) {
1076                ViewCell *viewcell = *it;
1077                if (!viewcell->Mailed()) {
1078                  viewcell->Mail();
1079                  currInfo.viewCellsBack++;
1080                }
1081                if (--viewcell->mCounter == 0)
[2053]1082                {
[492]1083                  currInfo.viewCellsFront--;
[2053]1084                }
[492]1085          }
[1199]1086#endif
1087         
[447]1088          break;
1089        }
1090        }
1091
1092        float position = (*ci).value;
1093       
1094        if (position > minBand && position < maxBand) {
[492]1095          currInfo.position = position;
[1199]1096
[492]1097         
[1199]1098          if (sumWeightsBack!=0.0f)
1099                info.contributionBack = sumContributionBack/sumWeightsBack;
1100          else
1101                info.contributionBack = 0.0f;
1102         
1103          if (sumWeightsFront!=0.0f)
1104                info.contributionFront = sumContributionFront/sumWeightsFront;
1105          else
1106                info.contributionFront = 0.0f;
1107         
[492]1108          GetCostRatio(
1109                                   leaf,
1110                                   currInfo);
1111         
[447]1112                       
1113          //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
1114          //      cout<<"cost= "<<sum<<endl;
[492]1115         
1116          if (currInfo.costRatio < info.costRatio) {
1117                info = currInfo;
[447]1118      }
1119    }
1120  }
1121 
[492]1122 
[447]1123  //  cout<<"===================="<<endl;
1124  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
1125  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
1126}
1127
1128void
[1233]1129RssTree::SortSubdivisionCandidates(
[447]1130                                                         RssTreeLeaf *node,
1131                                                         const int axis
1132                                                         )
1133{
1134 
1135  splitCandidates->clear();
1136 
[2130]1137  int requestedSize = 2*((int)node->rays.size());
[447]1138  // creates a sorted split candidates array
1139  if (splitCandidates->capacity() > 500000 &&
1140      requestedSize < (int)(splitCandidates->capacity()/10) ) {
1141   
1142    delete splitCandidates;
1143    splitCandidates = new vector<SortableEntry>;
1144  }
1145 
1146  splitCandidates->reserve(requestedSize);
1147
1148  // insert all queries
1149  for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
1150      ri < node->rays.end();
1151      ri++) {
1152        if ((*ri).mRay->IsActive()) {
1153          if (axis < 3) {
1154                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
[492]1155                                                                                                 (*ri).GetOrigin(axis),
1156                                                                                                 (void *)&(*ri))
[447]1157                                                                   );
1158          } else {
[492]1159                float pos = (*ri).GetDirParametrization(axis-3);
[447]1160                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
1161                                                                                                 pos,
[492]1162                                                                                                 (void *)&(*ri))
[447]1163                                                                   );
1164          }
1165        }
1166  }
1167 
1168  stable_sort(splitCandidates->begin(), splitCandidates->end());
1169}
1170
1171
1172void
1173RssTree::EvaluateLeafStats(const TraversalData &data)
1174{
1175
1176  // the node became a leaf -> evaluate stats for leafs
1177  RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
1178
1179  if (data.depth >= termMaxDepth)
1180    stat.maxDepthNodes++;
1181 
1182  //  if ( (int)(leaf->rays.size()) < termMinCost)
1183  //    stat.minCostNodes++;
[556]1184  if ( leaf->GetPvsSize() <= termMinPvs)
[447]1185        stat.minPvsNodes++;
1186
[556]1187  if ( leaf->GetPvsSize() <= termMinRays)
[447]1188        stat.minRaysNodes++;
1189
1190  if (leaf->GetAvgRayContribution() > termMaxRayContribution )
1191        stat.maxRayContribNodes++;
1192       
1193  if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
1194        stat.minSizeNodes++;
1195  }
1196
[2130]1197  if ((int)leaf->rays.size() > stat.maxRayRefs)
1198    stat.maxRayRefs = (int)leaf->rays.size();
[447]1199
1200}
1201
1202bool
1203RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
1204{
[556]1205  return ( (leaf->GetPvsSize() <= termMinPvs) ||
1206                   (leaf->rays.size() <= termMinRays) ||
[447]1207                   //                    (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
1208                   (leaf->depth >= termMaxDepth) ||
[464]1209                   (SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
[447]1210                   );
1211}
1212
1213
1214RssTreeNode *
1215RssTree::SubdivideNode(
1216                                           RssTreeLeaf *leaf,
1217                                           const AxisAlignedBox3 &box,
1218                                           AxisAlignedBox3 &backBBox,
1219                                           AxisAlignedBox3 &frontBBox
1220                                           )
1221{
1222 
1223  if (TerminationCriteriaSatisfied(leaf)) {
1224#if 0
1225        if (leaf->depth >= termMaxDepth) {
1226          cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
1227          cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
1228        }
1229#endif
1230               
1231        return leaf;
1232  }
1233       
[492]1234  SplitInfo info;
[447]1235       
1236  // select subdivision axis
[492]1237  int axis = SelectPlane( leaf,
1238                                                  box,
1239                                                  info
1240                                                  );
1241  //  Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
[447]1242 
1243  if (axis == -1) {
1244    return leaf;
1245  }
1246 
1247  stat.nodes+=2;
[569]1248  stat.leaves += 1;
[447]1249  stat.splits[axis]++;
1250
1251  // add the new nodes to the tree
1252  RssTreeInterior *node = new RssTreeInterior(leaf->parent);
1253
1254  node->axis = axis;
[492]1255  node->position = info.position;
[1867]1256
1257  node->bbox = GetBBox(leaf);
1258  node->dirBBox = GetDirBBox(leaf);
[447]1259 
[492]1260  RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
1261  RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);
[447]1262
[516]1263 
[492]1264 
[447]1265  // replace a link from node's parent
1266  if (  leaf->parent )
1267    leaf->parent->ReplaceChildLink(leaf, node);
1268  // and setup child links
1269  node->SetupChildLinks(back, front);
[516]1270
1271  back->bbox = leaf->bbox;
1272  front->bbox = leaf->bbox;
1273  back->dirBBox = leaf->dirBBox;
1274  front->dirBBox = leaf->dirBBox;
1275 
[447]1276  if (axis <= RssTreeNode::SPLIT_Z) {
[516]1277        back->bbox.SetMax(axis, info.position);
1278    front->bbox.SetMin(axis, info.position);
1279  } else {
1280        back->dirBBox.SetMax(axis-3, info.position);
1281    front->dirBBox.SetMin(axis-3, info.position);
[492]1282  }
[516]1283 
[492]1284  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1285          ri != leaf->rays.end();
1286          ri++) {
1287        if ((*ri).mRay->IsActive()) {
1288         
1289          // first unref ray from the former leaf
1290          (*ri).mRay->Unref();
1291         
1292          // Debug << "computed t: " << (*ri).mRay->mT << endl;
1293          // determine the side of this ray with respect to the plane
1294          int side = node->ComputeRaySide(*ri);
1295         
1296          if (side == 1)
1297                front->AddRay(*ri);
1298          else
1299                back->AddRay(*ri);
1300        } else
1301          (*ri).mRay->Unref();
[447]1302  }
1303
[492]1304  // distribute the total number of rays according to the distribution
1305  // of rays which remained
1306  //  front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
1307  //  back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
1308 
[447]1309#if 0
1310  front->SetPvsSize(pvsFront);
1311  back->SetPvsSize(pvsBack);
1312  // compute entropy as well
1313  front->ComputeEntropyImportance();
1314  back->ComputeEntropyImportance();
1315#else
[492]1316  UpdatePvsSize(front);
1317  UpdatePvsSize(back);
[447]1318#endif
1319
1320 
1321  // update stats
[485]1322  stat.rayRefs -= (int)leaf->rays.size();
[2130]1323  stat.rayRefs += (int)(back->rays.size() + front->rays.size());
[447]1324
1325 
1326  delete leaf;
1327  return node;
1328}
1329
1330
1331
1332
1333
1334
1335int
1336RssTree::ReleaseMemory(const int time)
1337{
1338  stack<RssTreeNode *> tstack;
1339 
1340  // find a node in the tree which subtree will be collapsed
1341  int maxAccessTime = time - accessTimeThreshold;
1342  int released;
1343
[516]1344  PushRoots(tstack);
[447]1345
1346  while (!tstack.empty()) {
1347    RssTreeNode *node = tstack.top();
1348    tstack.pop();
1349   
1350 
1351    if (!node->IsLeaf()) {
1352      RssTreeInterior *in = (RssTreeInterior *)node;
1353      //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1354      if (in->depth >= minCollapseDepth &&
1355                  in->lastAccessTime <= maxAccessTime) {
1356                released = CollapseSubtree(node, time);
1357                break;
1358      }
1359     
1360      if (in->back->GetAccessTime() < 
1361                  in->front->GetAccessTime()) {
1362                tstack.push(in->front);
1363                tstack.push(in->back);
1364      } else {
1365                tstack.push(in->back);
1366                tstack.push(in->front);
1367      }
1368    }
1369  }
1370
1371  while (tstack.empty()) {
1372    // could find node to collaps...
1373    //    cout<<"Could not find a node to release "<<endl;
1374    break;
1375  }
1376 
1377  return released;
1378}
1379
1380
1381
1382
1383RssTreeNode *
1384RssTree::SubdivideLeaf(
1385                                           RssTreeLeaf *leaf
1386                                           )
1387{
1388  RssTreeNode *node = leaf;
1389       
1390  AxisAlignedBox3 leafBBox = GetBBox(leaf);
1391
1392  static int pass = 0;
1393  pass ++;
1394       
1395  // check if we should perform a dynamic subdivision of the leaf
1396  if (!TerminationCriteriaSatisfied(leaf)) {
1397   
1398        // memory check and realese...
1399    if (GetMemUsage() > maxTotalMemory) {
1400      ReleaseMemory( pass );
1401    }
1402   
1403    AxisAlignedBox3 backBBox, frontBBox;
1404
1405    // subdivide the node
1406    node =
1407      SubdivideNode(leaf,
1408                                        leafBBox,
1409                                        backBBox,
1410                                        frontBBox
1411                                        );
1412  }
1413       
1414  return node;
1415}
1416
1417
1418
1419void
[492]1420RssTree::UpdateRays(
1421                                        VssRayContainer &remove,
[447]1422                                        VssRayContainer &add
1423                                        )
1424{
1425  RssTreeLeaf::NewMail();
1426
1427  // schedule rays for removal
1428  for(VssRayContainer::const_iterator ri = remove.begin();
1429      ri != remove.end();
1430      ri++) {
1431    (*ri)->ScheduleForRemoval();
1432  }
1433
1434  int inactive=0;
[1877]1435  Debug<<"Deleting rays..."<<endl<<flush;
[447]1436
1437  for(VssRayContainer::const_iterator ri = remove.begin();
1438      ri != remove.end();
1439      ri++) {
1440    if ((*ri)->ScheduledForRemoval())
1441          //      RemoveRay(*ri, NULL, false);
1442          // !!! BUG - with true it does not work correctly - aggreated delete
1443      RemoveRay(*ri, NULL, true);
1444    else
1445      inactive++;
1446  }
1447
[1877]1448  Debug<<"done."<<endl<<flush;
[447]1449
1450  //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
[1877]1451  Debug<<"Adding rays..."<<endl<<flush;
1452
[447]1453  for(VssRayContainer::const_iterator ri = add.begin();
1454      ri != add.end();
1455      ri++) {
1456        RssTreeNode::RayInfo info(*ri);
[516]1457        //      if (mForcedBoundingBox==NULL || ClipRay(info, bbox))
1458        AddRay(info);
[447]1459  }
[492]1460
[1877]1461  Debug<<"done."<<endl<<flush;
1462
[2130]1463  stat.rayRefs += (int)add.size() - (int)remove.size();
[492]1464
[1877]1465  Debug<<"Updating statistics..."<<endl<<flush;
[492]1466  UpdateTreeStatistics();
[1877]1467  Debug<<"done."<<endl<<flush;
1468
[492]1469  // check whether the tree should be prunned
1470  if (stat.rayRefs > mMaxRays) {
[1877]1471        Debug<<"Prunning rays..."<<endl<<flush;
[1906]1472        PruneRays((int)(0.8f*mMaxRays));
[1877]1473        Debug<<"done."<<endl<<flush;
[492]1474        //      UpdateTreeStatistics();
1475  }
[447]1476}
1477
[492]1478 
[447]1479
[492]1480
[447]1481void
1482RssTree::RemoveRay(VssRay *ray,
[2575]1483                   vector<RssTreeLeaf *> *affectedLeaves,
1484                   const bool removeAllScheduledRays
1485                   )
1486{       
1487  TstackEntry tstack;
[447]1488
[2575]1489  RssTreeNode::RayInfo ri(ray);
1490  PushRoots(tstack, ri);
[447]1491 
1492  RayTraversalData data;
1493
1494  // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1495
1496  while (!tstack.empty()) {
1497    data = tstack.top();
1498    tstack.pop();
1499
1500    if (!data.node->IsLeaf()) {
1501      // split the set of rays in two groups intersecting the
1502      // two subtrees
1503
1504      TraverseInternalNode(data, tstack);
1505     
1506    } else {
1507      // remove the ray from the leaf
1508      // find the ray in the leaf and swap it with the last ray...
1509      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
1510     
1511      if (!leaf->Mailed()) {
1512                leaf->Mail();
1513                if (affectedLeaves)
1514                  affectedLeaves->push_back(leaf);
1515       
1516                if (removeAllScheduledRays) {
[485]1517                  int tail = (int)leaf->rays.size()-1;
[447]1518
1519                  for (int i=0; i < (int)(leaf->rays.size()); i++) {
1520                        if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1521                          // find a ray to replace it with
1522                          while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1523                                stat.removedRayRefs++;
1524                                leaf->rays[tail].mRay->Unref();
1525                                leaf->rays.pop_back();
1526                                tail--;
1527                          }
1528
1529                          if (tail < i)
1530                                break;
1531             
1532                          stat.removedRayRefs++;
1533                          leaf->rays[i].mRay->Unref();
1534                          leaf->rays[i] = leaf->rays[tail];
1535                          leaf->rays.pop_back();
1536                          tail--;
1537                        }
1538                  }
1539                }
1540      }
1541     
1542      if (!removeAllScheduledRays)
1543                for (int i=0; i < (int)leaf->rays.size(); i++) {
1544                  if (leaf->rays[i].mRay == ray) {
1545                        stat.removedRayRefs++;
1546                        ray->Unref();
1547                        leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1548                        leaf->rays.pop_back();
1549                        // check this ray again
1550                        break;
1551                  }
1552                }
1553     
1554    }
1555  }
1556
1557  if (ray->RefCount() != 0) {
1558    cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1559    exit(1);
1560  }
1561 
1562}
1563
1564
1565void
1566RssTree::AddRay(RssTreeNode::RayInfo &info)
1567{
[1900]1568  if (EpsilonEqualV3(info.GetOrigin(), info.GetTermination(), Limits::Small))
1569        return;
[447]1570
1571  stack<RayTraversalData> tstack;
[516]1572
[569]1573  RssTreeNode *root = GetRoot(info.GetSourceObject());
[447]1574  tstack.push(RayTraversalData(root, info));
1575 
1576  RayTraversalData data;
1577 
[1900]1578
[447]1579  while (!tstack.empty()) {
1580    data = tstack.top();
1581    tstack.pop();
1582
1583    if (!data.node->IsLeaf()) {
1584      TraverseInternalNode(data, tstack);
1585    } else {
1586      // remove the ray from the leaf
1587      // find the ray in the leaf and swap it with the last ray...
1588      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
1589      leaf->AddRay(data.rayData);
1590      stat.addedRayRefs++;
1591    }
1592  }
1593}
1594
1595void
1596RssTree::TraverseInternalNode(
1597                                                          RayTraversalData &data,
1598                                                          stack<RayTraversalData> &tstack)
1599{
1600  RssTreeInterior *in =  (RssTreeInterior *) data.node;
1601 
[492]1602  // determine the side of this ray with respect to the plane
1603  int side = in->ComputeRaySide(data.rayData
1604                                                                );
[447]1605 
[492]1606  if (side == 1)
1607        tstack.push(RayTraversalData(in->front, data.rayData));
1608  else
1609        tstack.push(RayTraversalData(in->back, data.rayData));
1610 
[447]1611}
1612
1613
1614int
1615RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
1616{
1617  // first count all rays in the subtree
1618  // use mail 1 for this purpose
1619  stack<RssTreeNode *> tstack;
1620  int rayCount = 0;
1621  int totalRayCount = 0;
1622  int collapsedNodes = 0;
1623
1624#if DEBUG_COLLAPSE
1625  cout<<"Collapsing subtree"<<endl;
1626  cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1627  cout<<"depth="<<(int)sroot->depth<<endl;
1628#endif
1629
1630  //    tstat.collapsedSubtrees++;
1631  //    tstat.collapseDepths += (int)sroot->depth;
1632  //    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1633 
1634  tstack.push(sroot);
1635  VssRay::NewMail();
1636       
1637  while (!tstack.empty()) {
1638    collapsedNodes++;
1639    RssTreeNode *node = tstack.top();
1640    tstack.pop();
1641
1642    if (node->IsLeaf()) {
1643      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
1644      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1645                  ri != leaf->rays.end();
1646                  ri++) {
1647                               
1648                totalRayCount++;
1649                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
1650                  (*ri).mRay->Mail();
1651                  rayCount++;
1652                }
1653      }
1654    } else {
1655      tstack.push(((RssTreeInterior *)node)->back);
1656      tstack.push(((RssTreeInterior *)node)->front);
1657    }
1658  }
1659 
1660  VssRay::NewMail();
1661
1662  // create a new node that will hold the rays
1663  RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
1664  if (  newLeaf->parent )
1665    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1666 
1667
1668  tstack.push( sroot );
1669
1670  while (!tstack.empty()) {
1671
1672    RssTreeNode *node = tstack.top();
1673    tstack.pop();
1674
1675    if (node->IsLeaf()) {
1676      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
1677     
1678      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1679                  ri != leaf->rays.end();
1680                  ri++) {
1681                               
1682                // unref this ray from the old node
1683                               
1684                if ((*ri).mRay->IsActive()) {
1685                  (*ri).mRay->Unref();
1686                  if (!(*ri).mRay->Mailed()) {
1687                        (*ri).mRay->Mail();
1688                        newLeaf->AddRay(*ri);
1689                  }
1690                } else
1691                  (*ri).mRay->Unref();
1692                               
1693      }
1694    } else {
1695      tstack.push(((RssTreeInterior *)node)->back);
1696      tstack.push(((RssTreeInterior *)node)->front);
1697    }
1698  }
1699 
1700  // delete the node and all its children
1701  delete sroot;
1702 
1703  //   for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1704  //       ri != newLeaf->rays.end();
1705  //       ri++)
1706  //     (*ri).ray->UnMail(2);
1707
1708
1709#if DEBUG_COLLAPSE
1710  cout<<"Total memory before="<<GetMemUsage()<<endl;
1711#endif
1712
1713  stat.nodes -= collapsedNodes - 1;
[569]1714  stat.leaves -= collapsedNodes/2 - 1;
[447]1715  stat.rayRefs -= totalRayCount - rayCount;
1716 
1717#if DEBUG_COLLAPSE
1718  cout<<"collapsed nodes"<<collapsedNodes<<endl;
1719  cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1720  cout<<"Total memory after="<<GetMemUsage()<<endl;
1721  cout<<"================================"<<endl;
1722#endif
1723
1724  //  tstat.collapsedNodes += collapsedNodes;
1725  //  tstat.collapsedRays += totalRayCount - rayCount;
1726   
1727  return totalRayCount - rayCount;
1728}
1729
1730
1731int
1732RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1733{
1734  stack<RssTreeNode *> tstack;
[516]1735  PushRoots(tstack);
[447]1736
1737  Intersectable::NewMail();
1738  int pvsSize = 0;
1739       
1740  while (!tstack.empty()) {
1741    RssTreeNode *node = tstack.top();
1742    tstack.pop();
1743   
1744 
1745    if (node->IsLeaf()) {
1746          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1747          for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1748                  ri != leaf->rays.end();
1749                  ri++)
1750                if ((*ri).mRay->IsActive()) {
1751                  Intersectable *object;
[492]1752                  object = (*ri).GetObject();
[447]1753                  if (object && !object->Mailed()) {
1754                        pvsSize++;
1755                        object->Mail();
1756                  }
1757                }
1758        } else {
1759          RssTreeInterior *in = (RssTreeInterior *)node;
1760          if (in->axis < 3) {
1761                if (box.Max(in->axis) >= in->position )
1762                  tstack.push(in->front);
1763                               
1764                if (box.Min(in->axis) <= in->position )
1765                  tstack.push(in->back);
1766          } else {
1767                // both nodes for directional splits
1768                tstack.push(in->front);
1769                tstack.push(in->back);
1770          }
1771        }
1772  }
1773  return pvsSize;
1774}
1775
1776int
1777RssTree::CollectPvs(const AxisAlignedBox3 &box,
1778                                        ObjectContainer &pvs
1779                                        ) const
1780{
1781  stack<RssTreeNode *> tstack;
[516]1782  PushRoots(tstack);
[447]1783 
1784  Intersectable::NewMail();
1785 
1786  while (!tstack.empty()) {
1787    RssTreeNode *node = tstack.top();
1788    tstack.pop();
1789   
1790 
1791    if (node->IsLeaf()) {
1792          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1793          for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1794                  ri != leaf->rays.end();
1795                  ri++)
1796                if ((*ri).mRay->IsActive()) {
1797                  Intersectable *object;
[492]1798                  object = (*ri).GetObject();
[447]1799                  if (object && !object->Mailed()) {
1800                        pvs.push_back(object);
1801                        object->Mail();
1802                  }
1803                }
1804        } else {
1805          RssTreeInterior *in = (RssTreeInterior *)node;
1806          if (in->axis < 3) {
1807                if (box.Max(in->axis) >= in->position )
1808                  tstack.push(in->front);
1809               
1810                if (box.Min(in->axis) <= in->position )
1811                  tstack.push(in->back);
1812          } else {
1813                // both nodes for directional splits
1814                tstack.push(in->front);
1815                tstack.push(in->back);
1816          }
1817        }
1818  }
[2130]1819  return (int)pvs.size();
[447]1820}
1821
1822void
[492]1823RssTree::UpdateTreeStatistics()
[447]1824{
[1900]1825  for (int i=0; i < mRoots.size(); i++)
1826        UpdateImportance(mRoots[i]);
1827 
1828#if 0 
[447]1829  stack<RssTreeNode *> tstack;
1830       
1831  float sumPvsSize = 0.0f;
1832  float sumRayContribution = 0.0f;
[492]1833  float sumWeightedRayContribution = 0.0f;
[447]1834  float sumImportance = 0.0f;
1835  float sumPvsEntropy = 0.0f;
1836  float sumRayLengthEntropy = 0.0f;
[463]1837  float sumRays = 0.0f;
[447]1838 
1839  int leaves = 0;
1840
[1876]1841
1842  PushRoots(tstack);
[1900]1843 
[447]1844  while (!tstack.empty()) {
1845    RssTreeNode *node = tstack.top();
[1867]1846        tstack.pop();
1847       
[447]1848    if (node->IsLeaf()) {
1849          leaves++;
1850          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
[1867]1851          // updated above already
1852          //      UpdatePvsSize(leaf);
[463]1853         
[447]1854          sumPvsSize += leaf->GetPvsSize();
1855          sumRayContribution += leaf->GetAvgRayContribution();
[492]1856         
1857
1858          RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
1859          for (;it != leaf->rays.end(); ++it) {
1860                float weight, contribution;
1861                GetRayContribution(*it, weight, contribution);
1862                sumWeightedRayContribution += weight*contribution;
1863          }
1864         
1865          //      sumPvsEntropy += leaf->mPvsEntropy;
1866          //      sumRayLengthEntropy += leaf->mRayLengthEntropy;
[463]1867          sumRays += leaf->rays.size();
1868         
[447]1869          float imp = leaf->GetImportance();
1870         
[492]1871          //      if (imp > 1.0f)
1872          //            cout<<"warning imp > 1.0f:"<<imp<<endl;
[447]1873         
1874          sumImportance += imp;
1875         
1876        } else {
1877          RssTreeInterior *in = (RssTreeInterior *)node;
1878          tstack.push(in->front);
1879          tstack.push(in->back);
1880        }
1881  }
1882 
[492]1883  stat.avgPvsSize = sumPvsSize/(float)leaves;
1884  stat.avgRays = sumRays/(float)leaves;
1885  stat.avgRayContribution = sumRayContribution/(float)leaves;
1886  //  avgPvsEntropy = sumPvsEntropy/(float)leaves;
1887  //  avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
1888  stat.avgImportance = sumImportance/(float)leaves;
1889  stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
1890  stat.rayRefs = (int)sumRays;
[1900]1891#endif
[1867]1892 
[447]1893}
1894
[1867]1895void
1896RssTree::UpdateImportance(
1897                                                  RssTreeNode *node)
1898{
1899  if (node->IsLeaf()) {
1900        UpdatePvsSize((RssTreeLeaf *)node);
1901  } else {
1902        RssTreeInterior *in = (RssTreeInterior *)node;
1903        UpdateImportance(in->front);
1904        UpdateImportance(in->back);
1905        node->mImportance = in->front->mImportance + in->back->mImportance;
1906  }
1907}
[447]1908
[1867]1909// generate a single ray described by parameters in a unit cube
1910void
1911RssTree::GenerateRay(float *params, SimpleRayContainer &rays)
1912{
1913  RssTreeNode *node;
1914 
1915  // works only with single root now!
1916  node = mRoots[0];
1917  while (!node->IsLeaf()) {
1918        RssTreeInterior *in = (RssTreeInterior *)node;
1919        // decide whether to go left or right
1920        int axis = in->axis;
[1874]1921        //      float p = in->GetRelativePosition();
1922        float p = in->back->GetImportance()/(in->back->GetImportance() + in->front->GetImportance());
[1876]1923
1924        //      cout<<"p = "<<p<<endl;
[1867]1925        //      cout<<"axis="<<axis<<" pos="<<pos<<endl;
1926        //      cout<<GetBBox(in)<<endl;
1927        //      cout<<GetDirBBox(in)<<endl;
1928        //      cout<<"********"<<endl;
1929       
[1874]1930        if (params[axis] < p) {
[1867]1931          node = in->back;
[1876]1932          if (p > Limits::Small)
1933                params[axis] = params[axis] / p;
[1867]1934        } else {
1935          node = in->front;
[1876]1936          if (p < 1.0f - Limits::Small)
1937                params[axis] = (params[axis] - p) / (1.0f - p);
[1867]1938        }
1939  }
[492]1940
[1867]1941  //  cout<<params[0]<<" "<<params[1]<<" "<<params[2]<<" "<<params[3]<<" "<<params[4]<<endl;
1942 
1943  GenerateLeafRay(params, rays, node->bbox, node->dirBBox);
1944}
1945
[447]1946int
[492]1947RssTree::GenerateRays(
1948                                          const float ratioPerLeaf,
[447]1949                                          SimpleRayContainer &rays)
1950{
[1867]1951
[447]1952  stack<RssTreeNode *> tstack;
[516]1953  PushRoots(tstack);
[447]1954       
1955  while (!tstack.empty()) {
1956    RssTreeNode *node = tstack.top();
1957    tstack.pop();
1958               
1959    if (node->IsLeaf()) {
1960          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1961          float c = leaf->GetImportance();
[1189]1962          int num = int(c*ratioPerLeaf + 0.5f);
[447]1963          //                    cout<<num<<" ";
1964
1965          for (int i=0; i < num; i++) {
1966                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1967                Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1968                Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1969                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[1883]1970                rays.push_back(SimpleRay(origin,
1971                                                                 direction,
1972                                                                 SamplingStrategy::RSS_BASED_DISTRIBUTION,
1973                                                                 1.0f
1974                                                                 ));
[447]1975          }
1976        } else {
1977          RssTreeInterior *in = (RssTreeInterior *)node;
1978          // both nodes for directional splits
1979          tstack.push(in->front);
1980          tstack.push(in->back);
1981        }
1982  }
1983
[1563]1984  return (int)rays.size();
[447]1985}
1986
1987void
1988RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
1989{
1990  stack<RssTreeNode *> tstack;
[516]1991  PushRoots(tstack);
[447]1992       
1993  while (!tstack.empty()) {
1994    RssTreeNode *node = tstack.top();
1995    tstack.pop();
1996               
1997    if (node->IsLeaf()) {
1998          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1999          leaves.push_back(leaf);
2000        } else {
2001          RssTreeInterior *in = (RssTreeInterior *)node;
2002          // both nodes for directional splits
2003          tstack.push(in->front);
2004          tstack.push(in->back);
2005        }
2006  }
2007}
2008
2009bool
2010RssTree::ValidLeaf(RssTreeLeaf *leaf) const
2011{
[463]2012  return true;
2013  //return leaf->rays.size() > termMinRays/4;
[447]2014}
2015
2016
2017void
[464]2018RssTree::PickEdgeRays(RssTreeLeaf *leaf,
2019                                          int &r1,
2020                                          int &r2
2021                                          )
2022{
[2130]2023  int nrays = (int)leaf->rays.size();
[464]2024
2025  if (nrays == 2) {
2026        r1 = 0;
2027        r2 = 1;
2028        return;
2029  }
[492]2030
2031#if 1
[464]2032  int tries = min(20, nrays);
2033
2034  while (--tries) {
2035        r1 = Random(nrays);
2036        r2 = Random(nrays);
[492]2037        if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
[464]2038          break;
2039  }
2040
2041  if (r1 == r2)
[2130]2042        r2 = (r1+1)%(int)leaf->rays.size();
[492]2043
2044#else
2045  // pick a random ray
2046  int base = Random(nrays);
2047  RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
2048  Intersectable *baseObject = baseRay->GetObject();
2049 
2050  // and a random 5D derivative which will be used to find the minimal projected distances
2051  Vector3 spatialDerivative;
2052  Vector3 directionalDerivative;
2053
2054  while (1) {
2055        spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2056                                                                RandomValue(-1.0f, 1.0f),
2057                                                                RandomValue(-1.0f, 1.0f));
2058        float m = Magnitude(spatialDerivative);
2059        if (m != 0) {
2060          spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
2061          break;
2062        }
2063  }
2064
2065  while (1) {
2066        directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2067                                                                        RandomValue(-1.0f, 1.0f),
2068                                                                        0.0f);
2069        float m = Magnitude(directionalDerivative);
2070        if (m != 0) {
2071          directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
2072          break;
2073        }
2074  }
2075
2076  // now find the furthest sample from the same object and the closest from a different object
2077  int i = 0;
2078  float minDist = MAX_FLOAT;
2079  float maxDist = -MAX_FLOAT;
2080  r1 = base;
2081  r2 = base;
2082  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2083          ri != leaf->rays.end();
2084          ri++, i++) {
2085        float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
2086                                                                                                         *ri,
2087                                                                                                         spatialDerivative,
2088                                                                                                         directionalDerivative
2089                                                                                                         );
2090
2091        if ((*ri).GetObject() == baseObject) {
2092          if (dist > maxDist) {
2093                maxDist = dist;
2094                r1 = i;
2095          }
2096        } else {
2097          if (dist > 0.0f && dist < minDist) {
2098                minDist = dist;
2099                r2 = i;
2100          }
2101        }
2102  }
2103 
2104#endif
[464]2105}
2106
2107
[492]2108struct RayNeighbor {
2109  int rayInfo;
2110  float distance;
2111  RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
2112};
2113
[464]2114void
[492]2115RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
2116                                                        vector<RssTreeLeaf::SilhouetteRays> &rays
2117                                                        )
2118{
2119  // for every leaf find its neares neighbor from a different object
2120  vector<RayNeighbor> neighbors(leaf->rays.size());
2121  int i, j;
2122  for (i=0; i < leaf->rays.size(); i++)
2123        for (j=0; j < leaf->rays.size(); j++)
2124          if (i!=j) {
2125                float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
2126                if (d < neighbors[i].distance) {
2127                  neighbors[i].distance = d;
2128                  neighbors[i].rayInfo = j;
2129                }
2130          }
2131
2132  // now check which are the pairs of nearest neighbors
2133  for (i=0; i < leaf->rays.size(); i++) {
2134        int j = neighbors[i].rayInfo;
2135        if (neighbors[j].rayInfo == i) {
2136          // this is a silhouette edge pair
2137          if (i < j) {
2138                // generate an silhoutte ray pair
2139                rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
2140                                                                                                   &leaf->rays[j]));
2141          }
2142        } else {
2143          // this is not a silhouette ray - delete???
2144        }
2145  }
2146 
2147}
2148
2149
2150bool
2151GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
2152                                 const int index,
2153                                 int &i1,
2154                                 int &i2,
2155                                 int &i3)
2156{
2157  int found = 0;
2158  int indices[3];
2159
[2130]2160  int size = (int)rays.size();
[492]2161  // use russian roulete selection for the tripple
2162  // number of free positions for the bullet
2163  int positions = size - 1;
2164  int i;
2165  for (i=0; i < size; i++) {
2166        if (rays[i]->mRay->Mailed())
2167          positions--;
2168  }
2169 
2170  if (positions < 3)
2171        return false;
2172 
2173  for (i=0; i < size; i++) {
2174        if (i != index && !rays[i]->mRay->Mailed()) {
2175          float p = (3 - found)/(float)positions;
2176          if (Random(1.0f) < p) {
2177                indices[found] = i;
2178                found++;
2179          }
2180        }
2181        positions--;
2182  }
[500]2183  return true; // corr. matt
[492]2184}
2185
2186bool
2187RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
2188                                                                const RssTreeNode::RayInfo &r1,
2189                                                                const RssTreeNode::RayInfo &r2,
2190                                                                const RssTreeNode::RayInfo &r3)
2191{
2192 
2193
2194  return false;
2195}
2196
2197int
2198RssTree::RemoveInteriorRays(
2199                                                        RssTreeLeaf *leaf
2200                                                        )
2201{
2202#if 1
2203  // first collect all objects refered in this leaf
2204  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> > rayMap;
2205 
2206  RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
2207  for (; it != leaf->rays.end(); ++it) {
2208        Intersectable *object = (*it).GetObject();
2209       
2210        rayMap[object].push_back(&(*it));
2211//      vector<RayInfo *> *data = rayMap.Find(object);
2212//      if (data) {
2213//        data->push_back(&(*it));
2214//      } else {
2215//        //      rayMap[object] = vector<RayInfo *>;
2216//        rayMap[object].push_back(&(*it));
2217//      }
2218  }
2219
2220  // now go through all objects
2221  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> >::iterator mi;
2222
2223  // infos of mailed rays are scheduled for removal
2224  VssRay::NewMail();
2225  for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
2226        vector<RssTreeLeaf::RayInfo *> &rays = (*mi).second;
2227        vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
2228        int rayIndex = 0;
2229       
2230        for (; ri != rays.end(); ++ri, ++rayIndex) {
2231          RssTreeNode::RayInfo *ray = *ri;
[2130]2232          int tries = (int)rays.size();
[492]2233          for (int i = 0; i < tries; i++) {
2234                int r1, r2, r3;
2235                GetRandomTripple(rays,
2236                                                 rayIndex,
2237                                                 r1,
2238                                                 r2,
2239                                                 r3);
2240                if (IsRayConvexCombination(*ray,
2241                                                                   *rays[r1],
2242                                                                   *rays[r2],
2243                                                                   *rays[r3])) {
2244                  ray->mRay->Mail();
2245                }
2246          }
2247        }
2248  }
2249 
2250
2251#endif 
2252        return 0;
2253}
2254
[1715]2255
[492]2256void
[1715]2257RssTree::GenerateLeafRaysSimple(RssTreeLeaf *leaf,
2258                                                                const int numberOfRays,
2259                                                                SimpleRayContainer &rays)
2260{
2261 
2262  if (numberOfRays == 0)
2263        return;
2264 
2265  int startIndex = (int)rays.size();
[2119]2266   
[1715]2267  AxisAlignedBox3 box = GetBBox(leaf);
2268  AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2269
2270  float boxSize = Magnitude(box.Diagonal());
2271 
2272  const float smoothRange = 0.1f;
2273  if (smoothRange != 0.0f) {
2274        box.Scale(1.0f + smoothRange);
2275        dirBox.Scale(1.0f + smoothRange);
2276  }
2277
[1785]2278  //  Debug<<"R"<<flush;
[1715]2279 
2280  for (int i=0; i < numberOfRays; i++) {
2281        //      Debug<<i<<flush;
[1867]2282        float r[5];
[1715]2283       
[1867]2284        r[0] = RandomValue(0.0f, 1.0f);
2285        r[1] = RandomValue(0.0f, 1.0f);
2286        r[2] = RandomValue(0.0f, 1.0f);
2287        r[3] = RandomValue(0.0f, 1.0f);
2288        r[4] = RandomValue(0.0f, 1.0f);
[1715]2289       
[1867]2290        GenerateLeafRay(r, rays, box, dirBox);
[1715]2291  }
2292
[1785]2293  //  Debug<<"D"<<flush;
[1715]2294 
[2105]2295#if 0
[1715]2296  float density = 1.0f;
2297  float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
[1757]2298  for (int i=startIndex; i < rays.size(); i++) {
[1715]2299        rays[i].mPdf = p;
2300  }
[2105]2301#endif
[1715]2302}
2303
2304void
[1867]2305RssTree::GenerateLeafRay(
2306                                                 float *params,
2307                                                 SimpleRayContainer &rays,
[1876]2308                                                 const AxisAlignedBox3 &sbox,
2309                                                 const AxisAlignedBox3 &sDirBBox
[1867]2310                                                 )
2311{
[1876]2312  Vector3 origin = sbox.GetPoint(Vector3(params[0], params[1], params[2]));
2313  Vector3 dirVector = sDirBBox.GetPoint(Vector3(params[3], params[4], 0.0f));
[1867]2314  Vector3 direction = VssRay::GetInvDirParam(dirVector.x, dirVector.y);
2315 
2316  //    float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2317  float dist = 0.0f;
2318  float minT, maxT;
2319 
2320  // compute interection of the ray with the box
2321#if 1
[1876]2322  if (sbox.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
[1942]2323        dist = (maxT + Limits::Small);
[1867]2324#else
[1942]2325  float boxSize = Magnitude(sbox.Size());
[1867]2326  dist = 0.5f*boxSize;
2327#endif
2328 
2329  origin += direction*dist;
[1942]2330 
[1867]2331  //Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[1883]2332  rays.push_back(SimpleRay(origin,
2333                                                   direction,
2334                                                   SamplingStrategy::RSS_BASED_DISTRIBUTION,
2335                                                   1.0f
2336                                                   ));
[1867]2337}
2338
2339void
[447]2340RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
2341                                                  const int numberOfRays,
2342                                                  SimpleRayContainer &rays)
2343{
[463]2344
[492]2345  if (numberOfRays == 0)
2346        return;
2347 
[507]2348 
[1632]2349  int nrays = (int)leaf->rays.size();
2350  int startIndex = (int)rays.size();
[534]2351 
[492]2352  AxisAlignedBox3 box = GetBBox(leaf);
2353  AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2354
[507]2355
[1112]2356  const float smoothRange = 0.1f;
[534]2357  if (smoothRange != 0.0f) {
2358        box.Scale(1.0f + smoothRange);
2359        dirBox.Scale(1.0f + smoothRange);
[492]2360  }
2361 
2362  Exporter *exporter = NULL;
2363  VssRayContainer selectedRays;
2364  static int counter = 0;
2365  if (counter < 0) {
2366        char filename[128];
2367        sprintf(filename, "selected-rays%04d.x3d", counter);
2368        exporter = Exporter::GetExporter(filename);
2369        //      exporter->SetWireframe();
2370        //      exporter->ExportKdTree(*mKdTree);
2371        exporter->SetWireframe();
[1328]2372//      exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
[492]2373        exporter->SetWireframe();
2374        exporter->ExportBox(box);
2375        exporter->SetFilled();
2376        counter++;
2377  }
2378 
[1199]2379
2380  float extendedConvexCombinationProb = 0.0f; //0.7f
2381  //  float silhouetteCheckPercentage = 1.0f; //0.5f
2382  float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
2383  float silhouetteObjectCheckPercentage = 0.8f;
2384
[1112]2385  //  int numberOfTries = numberOfRays*4;
2386  int numberOfTries = numberOfRays;
[492]2387  int generated = 0;
2388
2389  vector<Intersectable *> pvs(0);
[1199]2390  if (silhouetteCheckPercentage != 0.0f) {
2391        Intersectable::NewMail();
2392        pvs.reserve(leaf->GetPvsSize());
2393        for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2394                ri != leaf->rays.end();
2395                ri++) {
2396          Intersectable *object = (*ri).GetObject();
2397          if (object) {
2398                if (!object->Mailed()) {
2399                  object->Mail();
[2053]2400                 
[1199]2401                  object->mCounter = 1;
2402                  pvs.push_back(object);
2403                  if (exporter)
2404                        exporter->ExportIntersectable(object);
2405                } else {
2406                  object->mCounter++;
2407                }
[556]2408          }
[492]2409        }
[1199]2410        // sort objects based on their mCounter
2411        sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
[492]2412  }
[556]2413 
[492]2414  for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
[549]2415        bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
2416                                                                                                                  extendedConvexCombinationProb));
[492]2417       
[447]2418        Vector3 origin, direction;
2419        // generate half of convex combination and half of random rays
2420        if (useExtendedConvexCombination) {
[492]2421          // pickup 2 random rays
[464]2422          int r1, r2;
2423          PickEdgeRays(leaf, r1, r2);
[492]2424
[463]2425         
2426          Vector3 o1 = leaf->rays[r1].GetOrigin();
2427          Vector3 o2 = leaf->rays[r2].GetOrigin();
2428         
[492]2429          const float overlap = 0.0f;
[464]2430          float w1, w2;
2431          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2432          origin = w1*o1 + w2*o2;
[492]2433          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
[447]2434          direction =
[492]2435                w1*leaf->rays[r1].GetDir() +
2436                w2*leaf->rays[r2].GetDir();
[1824]2437
2438          // shift the origin a little bit
2439          direction.Normalize();
2440         
[447]2441        } else {
[492]2442          Vector3 dirVector;
2443
2444          Vector3 pVector;
2445          Vector3 dVector;
2446
[1867]2447          pVector = Vector3(RandomValue(0.0f, 1.0f),
2448                                                RandomValue(0.0f, 1.0f),
2449                                                RandomValue(0.0f, 1.0f));
[492]2450         
[1867]2451          dVector = Vector3(RandomValue(0.0f, 1.0f),
2452                                                RandomValue(0.0f, 1.0f),
2453                                                0.0f);
[492]2454         
[1867]2455         
[534]2456          origin = box.GetPoint(pVector);
2457          dirVector = dirBox.GetPoint(dVector);
2458
[1824]2459          direction = VssRay::GetInvDirParam(dirVector.x,dirVector.y);
[447]2460        }
[492]2461       
[1824]2462       
[492]2463        //      float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2464        float dist = 0.0f;
2465        float minT, maxT;
2466        // compute interection of the ray with the box
2467
2468        if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2469          dist = maxT;
2470 
2471       
2472        origin += direction*dist;
2473       
2474        bool intersects = false;
2475       
[1112]2476        if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
[492]2477          if (exporter) {
2478                VssRay *ray = new VssRay(origin, origin + 100*direction, NULL, NULL);
2479                selectedRays.push_back(ray);
2480          }
[567]2481         
[492]2482          // check whether the ray does not intersect already visible objects
2483          Ray traversalRay;
2484          traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2485          for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2486                  oi != pvs.end();
2487                  oi++) {
2488                Intersectable *object = *oi;
[567]2489                // do not test every object
[1112]2490                if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
2491                         object->CastRay(traversalRay) ) {
[492]2492                  intersects = true;
2493                  break;
2494                }
2495          }
2496        }
2497       
2498       
2499        if (!intersects) {
2500          //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[1883]2501          rays.push_back(SimpleRay(origin,
2502                                                           direction,
2503                                                           SamplingStrategy::RSS_BASED_DISTRIBUTION,
2504                                                           1.0f
2505                                                           ));
[492]2506          generated++;
2507        }
[447]2508  }
[492]2509
[534]2510  //  cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
[556]2511  // assign pdfs to the generated rays
[1112]2512  // float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
[2105]2513#if 0
[1112]2514  float density = 1.0f;
[2105]2515 
[1112]2516  float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
[1757]2517  for (int i=startIndex; i < rays.size(); i++) {
[556]2518        rays[i].mPdf = p;
[537]2519  }
[2105]2520#endif
[537]2521 
[492]2522  if (exporter) {
2523        exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2524        delete exporter;
2525        CLEAR_CONTAINER(selectedRays);
2526  }
[534]2527 
[447]2528}
2529
2530int
[492]2531RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2532                                                 const float ratio)
2533{
2534  int i;
2535  int j;
2536 
2537  for (j=0, i=0; i < leaf->rays.size(); i++) {
2538        if (Random(1.0f) < ratio) {
2539          // copy a valid sample
2540          if (i!=j)
2541                leaf->rays[j] = leaf->rays[i];
2542          j++;
2543        } else {
2544          // delete the ray
2545          leaf->rays[i].mRay->Unref();
2546          if (leaf->rays[i].mRay->RefCount() != 0) {
2547                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2548                exit(1);
2549          }
2550          delete leaf->rays[i].mRay;
2551        }
2552  }
2553
2554
2555  leaf->rays.resize(j);
2556  int removed = (i-j);
2557  stat.rayRefs -= removed;
2558  return removed;
2559}
2560
2561int
2562RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2563                                                           const float ratio)
2564{
2565  int i;
2566 
2567  if (leaf->rays.size() == 0)
2568        return 0;
2569 
2570  sort(leaf->rays.begin(),
2571           leaf->rays.end(),
2572           RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2573
[1189]2574  int desired = int(ratio*leaf->rays.size());
2575  int removed = (int)leaf->rays.size() - desired;
[492]2576 
[1189]2577  for (i=desired; i < (int)leaf->rays.size(); i++) {
[492]2578        // delete the ray
2579        leaf->rays[i].mRay->Unref();
2580        if (leaf->rays[i].mRay->RefCount() != 0) {
2581          cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
[1632]2582          cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
[492]2583          exit(1);
2584        }
2585        delete leaf->rays[i].mRay;
2586  }
2587 
2588  leaf->rays.resize(desired);
2589  stat.rayRefs -= removed;
2590  return removed;
2591}
2592
2593
2594int
2595RssTree::PruneRays(
[1199]2596                                   const int maxRays
[492]2597                                   )
2598{
2599
2600  stack<RssTreeNode *> tstack;
2601  int prunned = 0;
2602
[1199]2603  // prune more rays to amortize  the procedure
[1501]2604  int desired = (int)(maxRays * 0.8f);
[492]2605
[1199]2606  Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
2607
[1877]2608  // prune random rays from each leaf so that the ratio's remain the same
2609  PushRoots(tstack);
2610 
2611  if ( maxRays >= stat.rayRefs )
2612        return 0;
2613 
2614  float ratio = desired/(float)stat.rayRefs;
2615 
2616  while (!tstack.empty()) {
2617        RssTreeNode *node = tstack.top();
2618        tstack.pop();
[492]2619       
[1877]2620        if (node->IsLeaf()) {
2621          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
[1883]2622          //prunned += PruneRaysRandom(leaf, ratio);
2623          prunned += PruneRaysContribution(leaf, ratio);
[1877]2624        } else {
2625          RssTreeInterior *in = (RssTreeInterior *)node;
2626          // both nodes for directional splits
2627          tstack.push(in->front);
2628          tstack.push(in->back);
[492]2629        }
2630  }
2631 
[1199]2632  Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
[492]2633 
2634  return prunned;
2635}
2636
2637int
[447]2638RssTree::GenerateRays(const int numberOfRays,
2639                                          const int numberOfLeaves,
2640                                          SimpleRayContainer &rays)
2641{
2642
[1931]2643  float param[5];
[1867]2644 
[1931]2645  for (int i=0; i < numberOfRays; i++) {
2646        mHalton.GetNext(5, param);
2647        GenerateRay(param, rays);
[1715]2648  }
2649 
[2130]2650  return (int)rays.size();
[1715]2651 
[447]2652}
2653
2654
2655float
2656RssTree::GetAvgPvsSize()
2657{
2658  stack<RssTreeNode *> tstack;
[516]2659  PushRoots(tstack);
[447]2660
2661  int sumPvs = 0;
2662  int leaves = 0;
2663  while (!tstack.empty()) {
2664    RssTreeNode *node = tstack.top();
2665    tstack.pop();
2666               
2667    if (node->IsLeaf()) {
2668          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2669          // update pvs size
[492]2670          UpdatePvsSize(leaf);
[447]2671          sumPvs += leaf->GetPvsSize();
2672          leaves++;
2673        } else {
2674          RssTreeInterior *in = (RssTreeInterior *)node;
2675          // both nodes for directional splits
2676          tstack.push(in->front);
2677          tstack.push(in->back);
2678        }
2679  }
2680
2681
2682  return sumPvs/(float)leaves;
2683}
2684
[1785]2685float weightAbsContributions = ABS_CONTRIBUTION_WEIGHT;
[492]2686// if small very high importance of the last sample
2687// if 1.0f then weighs = 1 1/2 1/3 1/4
[1891]2688float passSampleWeightDecay = 2.0f;
[1883]2689//float passSampleWeightDecay = 1000.0f;
[492]2690
[447]2691float
[492]2692RssTree::GetSampleWeight(const int pass)
[447]2693{
[492]2694  int passDiff = mCurrentPass - pass;
[1891]2695  if (passDiff < 0) {
2696        cerr<<"Error: passdif < 0"<<passDiff<<endl;
2697        exit(1);
2698  }
[492]2699  float weight;
[1891]2700  if (1)
[1876]2701        weight = 1.0f/(passDiff + passSampleWeightDecay);
[492]2702  else
[1876]2703        switch (passDiff) {
2704        case 0:
2705          //      weight = 1.0f;
2706          //      break;
2707        default:
2708          weight = 1.0f;
2709          break;
[492]2710          //    case 1:
[1876]2711          //      weight = 0.5f;
2712          //      break;
2713          //    case 2:
2714          //      weight = 0.25f;
2715          //      break;
2716          //    case 3:
2717          //      weight = 0.12f;
2718          //      break;
2719          //    case 4:
2720          //      weight = 0.06f;
2721          //      break;
2722          //    default:
2723          //      weight = 0.03f;
2724          //      break;
[492]2725        }
2726  return weight;
2727}
[447]2728
[492]2729void
[1112]2730RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2731                                                        float &weight,
2732                                                        float &contribution)
2733{
2734  VssRay *ray = info.mRay;
2735  weight = GetSampleWeight(mCurrentPass);
2736  contribution =
2737        weightAbsContributions*ray->mPvsContribution +
2738        (1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2739  // store the computed value
2740  info.mRay->mWeightedPvsContribution = weight*contribution;
2741}
2742
2743
2744void
[492]2745RssTree::ComputeImportance(RssTreeLeaf *leaf)
2746{
[1891]2747#if 0
[1877]2748  if (leaf->mTotalRays)
2749        leaf->mImportance = leaf->mTotalContribution / leaf->mTotalRays;
2750  else
[1883]2751        leaf->mImportance = Limits::Small;
[1877]2752  return;
[1883]2753#endif
[1877]2754 
[492]2755  if (0)
2756        leaf->mImportance = leaf->GetAvgRayContribution();
2757  else {
2758        RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2759          it_end = leaf->rays.end();
2760       
2761        float sumContributions = 0.0f;
2762        float sumRelContributions = 0.0f;
2763        float sumWeights = 0.0f;
[1824]2764        float maxContribution = 0.0f;
2765        float maxRelContribution = 0.0f;
2766        float sumLength = 0;
2767       
[492]2768        for (; it != it_end; ++it) {
2769          VssRay *ray = (*it).mRay;
[1824]2770          float sumLength = ray->Length();
2771
[492]2772          float weight = GetSampleWeight(ray->mPass);
2773         
[1824]2774          float c1 = weight*ray->mPvsContribution;
2775          float c2 = weight*ray->mRelativePvsContribution;
[1900]2776
[1824]2777          sumContributions += c1;
2778          if (c1 > maxContribution)
[1942]2779                maxContribution = c1;
[1824]2780         
2781          //$$ 20.7. changed to sqr to pronouce higher contribution so that
2782          // they do not get smoothed!!
[1199]2783          //sumRelContributions += weight*ray->mRelativePvsContribution;
[1824]2784         
2785          sumRelContributions += c2;
2786          if (c2 > maxRelContribution)
[1900]2787                maxRelContribution = c2;
2788         
[1112]2789          //sumWeights += weight;
[492]2790          sumWeights += 1.0f;
2791        }
[1824]2792
2793        float distImportance = 0.0f;
2794
[1829]2795        float spatialArea = GetBBox(leaf).SurfaceArea();
2796        float dirArea = GetDirBBox(leaf).SurfaceArea();
[1824]2797        if (leaf->rays.size()) {
2798          // compute average length of the ray
2799          float avgLength = sumLength/leaf->rays.size();
2800          // compute estimated area of the visible surface corresponding to these rays
[1829]2801          float ss = spatialArea/2.0f;
2802          float sd = dirArea;
[1824]2803          float s = ss + avgLength*(2*sqrt(ss*sd) + avgLength*sd);
2804         
2805          distImportance = s/leaf->rays.size();
2806        }
2807       
[492]2808        // $$
[1867]2809        // sumWeights = leaf->mTotalRays;
[1824]2810
[1829]2811        //$$ test 30.11. 2006
2812        // normalize per volume of the cell not the number of rays!
[1832]2813        // sumWeights = spatialArea*dirArea;
[1829]2814       
[1824]2815        if (sumWeights != 0.0f) {
[1891]2816#if 1
[1824]2817          leaf->mImportance =
[1891]2818                pow(((weightAbsContributions*sumContributions +
[1900]2819                          (1.0f - weightAbsContributions)*sumRelContributions)/sumWeights), 1.0f);
[1942]2820#else
[1891]2821          leaf->mImportance =
2822                (weightAbsContributions*maxContribution +
2823                 (1.0f - weightAbsContributions)*maxRelContribution)/sumWeights;
2824#endif
[1824]2825        } else
[492]2826          leaf->mImportance = 0.0f;
2827       
[447]2828  }
[492]2829 
2830  // return GetAvgRayContribution()*mImportance;
2831  //return GetAvgRayContribution();
[447]2832}
2833
2834
2835float
[1867]2836RssTreeNode::GetImportance()  const
[447]2837{
[1891]2838  // return mImportance*mImportance;
[492]2839  return mImportance;
2840}
2841
2842
2843float
2844RssTreeLeaf::ComputePvsEntropy() const
2845{
[447]2846  int samples = 0;
2847  Intersectable::NewMail();
2848  // set all object as belonging to the fron pvs
[492]2849  for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
[447]2850          ri != rays.end();
2851          ri++)
2852        if ((*ri).mRay->IsActive()) {
[492]2853          Intersectable *object = (*ri).GetObject();
[447]2854          if (object) {
2855                if (!object->Mailed()) {
2856                  object->Mail();
2857                  object->mCounter = 1;
2858                } else
[2053]2859                {
2860               
[447]2861                  object->mCounter++;
[2053]2862                }
[447]2863                samples++;
2864          }
2865        }
2866 
2867  float entropy = 0.0f;
2868 
2869  if (samples > 1) {
2870        Intersectable::NewMail();
2871        for(RayInfoContainer::const_iterator ri = rays.begin();
2872                ri != rays.end();
2873                ri++)
2874          if ((*ri).mRay->IsActive()) {
[492]2875                Intersectable *object = (*ri).GetObject();
[447]2876                if (object) {
2877                  if (!object->Mailed()) {
2878                        object->Mail();
2879                        float p = object->mCounter/(float)samples;
2880                        entropy -= p*log(p);
2881                  }
2882                }
2883          }
2884        entropy = entropy/log((float)samples);
2885  }
2886  else
2887        entropy = 1.0f;
2888 
2889  return entropy;
2890}
2891
2892float
[492]2893RssTreeLeaf::ComputeRayLengthEntropy() const
[447]2894{
2895  // get sum of all ray lengths
2896  // consider only passing rays or originating rays
2897  float sum = 0.0f;
2898  int samples = 0;
2899  int i=0;
2900  for(RayInfoContainer::const_iterator ri = rays.begin();
2901      ri != rays.end();
2902      ri++)
2903        if ((*ri).mRay->IsActive()) {
2904//        float s;
2905//        if (i == 0)
2906//              s = 200;
2907//        else
2908//              s = 100;
2909//        i++;
2910         
2911          sum += (*ri).mRay->GetSize();
2912          samples++;
2913        }
2914 
2915 
2916  float entropy = 0.0f;
2917 
2918  if (samples > 1) {
2919        i = 0;
2920        for(RayInfoContainer::const_iterator ri = rays.begin();
2921                ri != rays.end();
2922                ri++)
2923          if ((*ri).mRay->IsActive()) {
2924//              float s;
2925//              if (i==0)
2926//                s = 200;
2927//              else
2928//                s = 100;
2929//              i++;
2930//              float p = s/sum;
2931                float p = (*ri).mRay->GetSize()/sum;
2932                entropy -= p*log(p);
2933          }
2934        entropy = entropy/log((float)samples);
2935  } else
2936        entropy = 1.0f;
2937 
2938  return entropy;
2939}
2940 
2941
2942
[492]2943float
2944RssTreeLeaf::ComputeEntropyImportance() const
[447]2945{
2946 
2947  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
[492]2948  return 1.0f - ComputePvsEntropy();
2949}
2950
2951float
2952RssTreeLeaf::ComputeRayContributionImportance() const
2953{
2954  //  mPvsEntropy = ComputePvsEntropy();
2955  //  mRayLengthEntropy = ComputeRayLengthEntropy();
[447]2956 
[492]2957  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
2958  return 1.0f - ComputePvsEntropy();
[447]2959}
2960
2961
2962AxisAlignedBox3
2963RssTree::GetShrankedBBox(const RssTreeNode *node) const
2964{
2965  if (node->parent == NULL)
2966        return bbox;
2967 
2968  if (!node->IsLeaf())
2969        return ((RssTreeInterior *)node)->bbox;
2970
2971  // evaluate bounding box from the ray origins
2972  AxisAlignedBox3 box;
2973  box.Initialize();
2974  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2975  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2976          ri != leaf->rays.end();
2977          ri++)
2978        if ((*ri).mRay->IsActive()) {
2979          box.Include((*ri).GetOrigin());
2980        }
2981  return box;
2982}
[466]2983
2984
2985int
2986RssTree::CollectRays(VssRayContainer &rays,
2987                                         const int number)
2988{
2989  VssRayContainer allRays;
2990  CollectRays(allRays);
[492]2991
[466]2992 
2993  int desired = min(number, (int)allRays.size());
2994  float prob = desired/(float)allRays.size();
[1581]2995  //int skip = allRays.size()/desired;
[466]2996  while (rays.size() < desired) {
2997        VssRayContainer::const_iterator it = allRays.begin();
2998        for (; it != allRays.end() && rays.size() < desired; it++) {
2999          if (Random(1.0f) < prob)
3000                rays.push_back(*it);
3001        }
3002  }
[485]3003  return (int)rays.size();
[466]3004}
3005
3006
3007int
3008RssTree::CollectRays(VssRayContainer &rays
3009                                         )
3010{
3011  VssRay::NewMail();
3012
3013  stack<RssTreeNode *> tstack;
[516]3014  PushRoots(tstack);
[466]3015 
3016  while (!tstack.empty()) {
3017    RssTreeNode *node = tstack.top();
3018    tstack.pop();
3019       
3020    if (node->IsLeaf()) {
3021          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
3022          // update pvs size
3023          RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
3024          for (;it != leaf->rays.end(); ++it)
3025                if (!(*it).mRay->Mailed()) {
3026                  (*it).mRay->Mail();
3027                  rays.push_back((*it).mRay);
3028                }
3029        } else {
3030          RssTreeInterior *in = (RssTreeInterior *)node;
3031          // both nodes for directional splits
3032          tstack.push(in->front);
3033          tstack.push(in->back);
3034        }
3035  }
3036 
[485]3037  return (int)rays.size();
[466]3038}
[492]3039
[516]3040
3041RssTreeNode *
3042RssTree::GetRoot(Intersectable *object) const
3043{
3044  if (mPerObjectTree && object) {
[569]3045        int id = object->GetId()-1;
3046        if (id < 0 || id >= mRoots.size()) {
[1632]3047          Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
[569]3048                endl<<flush;
[579]3049          id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
[569]3050        }
[516]3051        return mRoots[id];
3052  } else
3053        return mRoots[0];
3054}
3055
3056void
3057RssTree::PushRoots(priority_queue<TraversalData> &st) const
3058{
3059  for (int i=0; i < mRoots.size(); i++)
3060        st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
3061}
3062
3063void
3064RssTree::PushRoots(stack<RssTreeNode *> &st)  const
3065{
3066  for (int i=0; i < mRoots.size(); i++)
3067        st.push(mRoots[i]);
3068}
3069
3070void
3071RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info)  const
3072{
3073  for (int i=0; i < mRoots.size(); i++)
3074        st.push(RayTraversalData(mRoots[i], info));
3075}
[860]3076
3077
[1867]3078
[1884]3079bool
3080RssBasedDistribution::GenerateSample(SimpleRay &ray)
3081{
3082  float r[5];
[1867]3083
[1884]3084  mHalton.GetNext(5, r);
[1867]3085
[1891]3086  if (mRssTree == NULL) {
3087        // use direction based distribution
3088        Vector3 origin, direction;
3089        mPreprocessor.mViewCellsManager->GetViewPoint(origin, Vector3(r[0], r[1], r[2]));
3090       
3091        direction = UniformRandomVector(r[3], r[4]);
3092
3093        const float pdf = 1.0f;
3094        ray = SimpleRay(origin, direction, RSS_BASED_DISTRIBUTION, pdf);
3095       
3096  } else {
3097       
3098        SimpleRayContainer rays;
3099        mRssTree->GenerateRay(r, rays);
3100        ray = rays[0];
3101  }
[1884]3102 
[1891]3103  return true;
[1112]3104}
[1884]3105
3106void
3107RssBasedDistribution::Update(VssRayContainer &vssRays)
3108{
3109  if (mRssTree == NULL) {
3110        mRssTree = new RssTree;
3111        mRssTree->SetPass(mPass);
3112       
3113        /// compute view cell contribution of rays if view cells manager already constructed
3114        //  mViewCellsManager->ComputeSampleContributions(mVssRays, true, false);
3115        mRssTree->Construct(mPreprocessor.mObjects, vssRays);
3116        mRssTree->stat.Print(mPreprocessor.mStats);
3117        cout<<"RssTree root PVS size = "<<mRssTree->GetRootPvsSize()<<endl;
3118       
3119  } else {
[1891]3120        mRssTree->SetPass(mPass);
[1884]3121        Debug<<"Adding rays...\n"<<flush;
3122        mRssTree->AddRays(vssRays);
3123        Debug<<"done.\n"<<flush;
3124        if (1) {
3125          //    if (mUpdateSubdivision) {
3126          Debug<<"Updating rss tree...\n"<<flush;
3127          int subdivided = mRssTree->UpdateSubdivision();
3128          Debug<<"done.\n"<<flush;
3129          cout<<"subdivided leafs = "<<subdivided<<endl;
3130          cout<<"#total leaves = "<<mRssTree->stat.Leaves()<<endl;
3131        }
3132  }
3133  mPass++;
3134}
3135
3136}
Note: See TracBrowser for help on using the repository browser.