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

Revision 1824, 73.7 KB checked in by bittner, 18 years ago (diff)

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