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

Revision 1112, 66.7 KB checked in by bittner, 18 years ago (diff)

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