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

Revision 863, 65.5 KB checked in by mattausch, 19 years ago (diff)

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