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

Revision 1942, 74.3 KB checked in by bittner, 18 years ago (diff)

tmp commit

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