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

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