source: trunk/VUT/GtpVisibilityPreprocessor/src/RssTree.cpp @ 563

Revision 563, 65.2 KB checked in by bittner, 18 years ago (diff)

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