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

Revision 537, 64.8 KB checked in by bittner, 19 years ago (diff)

ray probability density support
Line 
1// ================================================================
2// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
3// ****************************************************************
4/**
5   The KD tree based LSDS
6*/
7// Initial coding by
8/**
9   @author Jiri Bittner
10*/
11
12// Standard headers
13#include <stack>
14#include <queue>
15#include <algorithm>
16#include <fstream>
17#include <string>
18
19#include "RssTree.h"
20
21#include "Environment.h"
22#include "VssRay.h"
23#include "Intersectable.h"
24#include "Ray.h"
25#include "Containers.h"
26#include "ViewCell.h"
27#include "Exporter.h"
28#include "Preprocessor.h"
29#include "SceneGraph.h"
30
31#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 minCost )\n"<<
222    minPvsNodes*100/(double)Leaves()<<endl;
223
224  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\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 && !object->Mailed()) {
2084          object->Mail();
2085          pvs.push_back(object);
2086          if (exporter)
2087                exporter->ExportIntersectable(object);
2088        }
2089  }
2090 
2091  for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2092        bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) < 0.7f));
2093       
2094        Vector3 origin, direction;
2095        // generate half of convex combination and half of random rays
2096        if (useExtendedConvexCombination) {
2097          // pickup 2 random rays
2098          int r1, r2;
2099          PickEdgeRays(leaf, r1, r2);
2100
2101         
2102          Vector3 o1 = leaf->rays[r1].GetOrigin();
2103          Vector3 o2 = leaf->rays[r2].GetOrigin();
2104         
2105          const float overlap = 0.0f;
2106          float w1, w2;
2107          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2108          origin = w1*o1 + w2*o2;
2109          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2110          direction =
2111                w1*leaf->rays[r1].GetDir() +
2112                w2*leaf->rays[r2].GetDir();
2113        } else {
2114          Vector3 dirVector;
2115
2116          Vector3 pVector;
2117          Vector3 dVector;
2118
2119          bool useHalton = true;
2120         
2121          if (useHalton) {
2122                // generate a random 5D vector
2123                pVector = Vector3(leaf->halton.GetNumber(1),
2124                                                  leaf->halton.GetNumber(2),
2125                                                  leaf->halton.GetNumber(3));
2126               
2127                dVector = Vector3(leaf->halton.GetNumber(4),
2128                                                  leaf->halton.GetNumber(5),
2129                                                  0.0f);
2130               
2131                leaf->halton.GenerateNext();
2132          } else {
2133                pVector = Vector3(RandomValue(0.0f, 1.0f),
2134                                                  RandomValue(0.0f, 1.0f),
2135                                                  RandomValue(0.0f, 1.0f));
2136               
2137                dVector = Vector3(RandomValue(0.0f, 1.0f),
2138                                                  RandomValue(0.0f, 1.0f),
2139                                                  0.0f);
2140          }
2141         
2142          origin = box.GetPoint(pVector);
2143          dirVector = dirBox.GetPoint(dVector);
2144
2145          direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
2146        }
2147
2148        // shift the origin a little bit
2149        direction.Normalize();
2150       
2151        //      float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2152        float dist = 0.0f;
2153        float minT, maxT;
2154        // compute interection of the ray with the box
2155
2156        if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2157          dist = maxT;
2158 
2159       
2160        origin += direction*dist;
2161       
2162        bool intersects = false;
2163       
2164
2165        if (i > numberOfRays/2) {
2166          if (exporter) {
2167                VssRay *ray = new VssRay(origin, origin + 100*direction, NULL, NULL);
2168                selectedRays.push_back(ray);
2169          }
2170
2171          // check whether the ray does not intersect already visible objects
2172          Ray traversalRay;
2173          traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2174          for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2175                  oi != pvs.end();
2176                  oi++) {
2177                Intersectable *object = *oi;
2178               
2179                if ( object->CastRay(traversalRay) ) {
2180                  intersects = true;
2181                  break;
2182                }
2183          }
2184        }
2185       
2186       
2187        if (!intersects) {
2188          //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2189          rays.push_back(SimpleRay(origin, direction));
2190          generated++;
2191        }
2192  }
2193
2194 
2195  //  cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2196  // assign probabilitites to the generated rays
2197  float p = (box.GetVolume()*dirBox.GetVolume())/generated;
2198  for (i=startIndex; i < rays.size(); i++) {
2199        rays[i].mProbability = p;
2200        rays.mSumProbabilities += p;
2201  }
2202 
2203 
2204  if (exporter) {
2205        exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2206        delete exporter;
2207        CLEAR_CONTAINER(selectedRays);
2208  }
2209 
2210 
2211}
2212
2213int
2214RssTree::PruneRays(RssTreeLeaf *leaf,
2215                                   const float contributionThreshold)
2216{
2217  int i;
2218  int j;
2219 
2220  for (j=0, i=0; i < leaf->rays.size(); i++) {
2221       
2222        if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
2223          // copy a valid sample
2224          if (i!=j)
2225                leaf->rays[j] = leaf->rays[i];
2226          j++;
2227        } else {
2228          // delete the ray
2229          leaf->rays[i].mRay->Unref();
2230          if (leaf->rays[i].mRay->RefCount() != 0) {
2231                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2232                exit(1);
2233          }
2234          delete leaf->rays[i].mRay;
2235        }
2236  }
2237
2238
2239  leaf->rays.resize(j);
2240  int removed = (i-j);
2241  stat.rayRefs -= removed;
2242  return removed;
2243}
2244
2245int
2246RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2247                                                 const float ratio)
2248{
2249  int i;
2250  int j;
2251 
2252  for (j=0, i=0; i < leaf->rays.size(); i++) {
2253       
2254        if (Random(1.0f) < ratio) {
2255          // copy a valid sample
2256          if (i!=j)
2257                leaf->rays[j] = leaf->rays[i];
2258          j++;
2259        } else {
2260          // delete the ray
2261          leaf->rays[i].mRay->Unref();
2262          if (leaf->rays[i].mRay->RefCount() != 0) {
2263                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2264                exit(1);
2265          }
2266          delete leaf->rays[i].mRay;
2267        }
2268  }
2269
2270
2271  leaf->rays.resize(j);
2272  int removed = (i-j);
2273  stat.rayRefs -= removed;
2274  return removed;
2275}
2276
2277int
2278RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2279                                                           const float ratio)
2280{
2281  int i;
2282 
2283  if (leaf->rays.size() == 0)
2284        return 0;
2285 
2286  sort(leaf->rays.begin(),
2287           leaf->rays.end(),
2288           RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2289
2290  int desired = ratio*leaf->rays.size();
2291  int removed = leaf->rays.size() - desired;
2292 
2293  for (i=desired; i < leaf->rays.size(); i++) {
2294        // delete the ray
2295        leaf->rays[i].mRay->Unref();
2296        if (leaf->rays[i].mRay->RefCount() != 0) {
2297          cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2298          exit(1);
2299        }
2300        delete leaf->rays[i].mRay;
2301  }
2302 
2303  leaf->rays.resize(desired);
2304  stat.rayRefs -= removed;
2305  return removed;
2306}
2307
2308
2309int
2310RssTree::PruneRays(
2311                                   const int desired
2312                                   )
2313{
2314  bool globalPrunning = false;
2315
2316  stack<RssTreeNode *> tstack;
2317  int prunned = 0;
2318
2319  Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl;
2320
2321  if (globalPrunning) {
2322        VssRayContainer allRays;
2323        allRays.reserve(stat.rayRefs);
2324        CollectRays(allRays);
2325        sort(allRays.begin(),
2326                 allRays.end(),
2327                 GreaterWeightedPvsContribution);
2328       
2329        if ( desired >= allRays.size() )
2330          return 0;
2331       
2332        float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
2333       
2334        PushRoots(tstack);
2335       
2336        while (!tstack.empty()) {
2337          RssTreeNode *node = tstack.top();
2338          tstack.pop();
2339         
2340          if (node->IsLeaf()) {
2341                RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2342                prunned += PruneRays(leaf, contributionThreshold);
2343          } else {
2344                RssTreeInterior *in = (RssTreeInterior *)node;
2345                // both nodes for directional splits
2346                tstack.push(in->front);
2347                tstack.push(in->back);
2348          }
2349        }
2350  } else {
2351        // prune random rays from each leaf so that the ratio's remain the same
2352        PushRoots(tstack);
2353        float ratio = desired/(float)stat.rayRefs;
2354       
2355        while (!tstack.empty()) {
2356          RssTreeNode *node = tstack.top();
2357          tstack.pop();
2358         
2359          if (node->IsLeaf()) {
2360                RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2361                //              prunned += PruneRaysRandom(leaf, ratio);
2362                prunned += PruneRaysContribution(leaf, ratio);
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
2371       
2372
2373
2374  }
2375 
2376 
2377 
2378  Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl;
2379 
2380  return prunned;
2381}
2382
2383int
2384RssTree::GenerateRays(const int numberOfRays,
2385                                          const int numberOfLeaves,
2386                                          SimpleRayContainer &rays)
2387{
2388
2389
2390  vector<RssTreeLeaf *> leaves;
2391       
2392  CollectLeaves(leaves);
2393       
2394  sort(leaves.begin(),
2395           leaves.end(),
2396           GreaterContribution);
2397                         
2398
2399  float sumContrib = 0.0;
2400  int i;
2401  int k = 0;
2402  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2403        if (ValidLeaf(leaves[i])) {
2404          float c = leaves[i]->GetImportance();
2405          sumContrib += c;
2406          //            cout<<"ray contrib "<<i<<" : "<<c<<endl;
2407          k++;
2408        }
2409                               
2410  float avgContrib = sumContrib/numberOfLeaves;
2411 
2412  // always generate at leat n ray per leaf
2413  int fixedPerLeaf = 1;
2414  int fixed = 1*leaves.size();
2415  int iGenerated = numberOfRays;
2416  float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2417
2418  k = 0;
2419  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2420        if (ValidLeaf(leaves[i])) {
2421          k++;
2422          RssTreeLeaf *leaf = leaves[i];
2423          float c = leaf->GetImportance();
2424
2425          int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2426          GenerateLeafRays(leaf, num, rays);
2427        }
2428
2429  return rays.size();
2430}
2431
2432
2433float
2434RssTree::GetAvgPvsSize()
2435{
2436  stack<RssTreeNode *> tstack;
2437  PushRoots(tstack);
2438
2439  int sumPvs = 0;
2440  int leaves = 0;
2441  while (!tstack.empty()) {
2442    RssTreeNode *node = tstack.top();
2443    tstack.pop();
2444               
2445    if (node->IsLeaf()) {
2446          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2447          // update pvs size
2448          UpdatePvsSize(leaf);
2449          sumPvs += leaf->GetPvsSize();
2450          leaves++;
2451        } else {
2452          RssTreeInterior *in = (RssTreeInterior *)node;
2453          // both nodes for directional splits
2454          tstack.push(in->front);
2455          tstack.push(in->back);
2456        }
2457  }
2458
2459
2460  return sumPvs/(float)leaves;
2461}
2462
2463float weightAbsContributions = 0.0f;
2464// if small very high importance of the last sample
2465// if 1.0f then weighs = 1 1/2 1/3 1/4
2466float passSampleWeightDecay = 0.5f;
2467
2468void
2469RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2470                                                        float &weight,
2471                                                        float &contribution)
2472{
2473  VssRay *ray = info.mRay;
2474  weight = 1.0f/(mCurrentPass - ray->mPass + passSampleWeightDecay);
2475  contribution =
2476        weightAbsContributions*ray->mPvsContribution +
2477        (1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2478  // store the computed value
2479  info.mRay->mWeightedPvsContribution = weight*contribution;
2480}
2481
2482float
2483RssTree::GetSampleWeight(const int pass)
2484{
2485  int passDiff = mCurrentPass - pass;
2486  float weight;
2487  if (1)
2488        weight = 1.0f/(passDiff + passSampleWeightDecay);
2489  else
2490        switch (passDiff) {
2491        case 0:
2492          weight = 1.0f;
2493          break;
2494        default:
2495          weight = 0.01f;
2496          break;
2497          //    case 1:
2498//        weight = 0.5f;
2499//        break;
2500//      case 2:
2501//        weight = 0.25f;
2502//        break;
2503//      case 3:
2504//        weight = 0.12f;
2505//        break;
2506//      case 4:
2507//        weight = 0.06f;
2508//        break;
2509//      default:
2510//        weight = 0.03f;
2511//        break;
2512        }
2513  return weight;
2514}
2515
2516void
2517RssTree::ComputeImportance(RssTreeLeaf *leaf)
2518{
2519  if (0)
2520        leaf->mImportance = leaf->GetAvgRayContribution();
2521  else {
2522        RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2523          it_end = leaf->rays.end();
2524       
2525        float sumContributions = 0.0f;
2526        float sumRelContributions = 0.0f;
2527        float sumWeights = 0.0f;
2528        for (; it != it_end; ++it) {
2529          VssRay *ray = (*it).mRay;
2530          float weight = GetSampleWeight(ray->mPass);
2531         
2532          sumContributions += weight*ray->mPvsContribution;
2533          sumRelContributions += weight*ray->mRelativePvsContribution;
2534         
2535          //      sumWeights += weight;
2536          sumWeights += 1.0f;
2537        }
2538        // $$
2539        // sumWeights = leaf->mTotalRays;
2540         
2541        if (sumWeights != 0.0f)
2542          leaf->mImportance =
2543                (weightAbsContributions*sumContributions +
2544                 (1.0f - weightAbsContributions)*sumRelContributions)/sumWeights;
2545        else
2546          leaf->mImportance = 0.0f;
2547       
2548  }
2549 
2550  // return GetAvgRayContribution()*mImportance;
2551  //return GetAvgRayContribution();
2552}
2553
2554
2555float
2556RssTreeLeaf::GetImportance()  const
2557{
2558  return mImportance;
2559}
2560
2561
2562float
2563RssTreeLeaf::ComputePvsEntropy() const
2564{
2565  int samples = 0;
2566  Intersectable::NewMail();
2567  // set all object as belonging to the fron pvs
2568  for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2569          ri != rays.end();
2570          ri++)
2571        if ((*ri).mRay->IsActive()) {
2572          Intersectable *object = (*ri).GetObject();
2573          if (object) {
2574                if (!object->Mailed()) {
2575                  object->Mail();
2576                  object->mCounter = 1;
2577                } else
2578                  object->mCounter++;
2579                samples++;
2580          }
2581        }
2582 
2583  float entropy = 0.0f;
2584 
2585  if (samples > 1) {
2586        Intersectable::NewMail();
2587        for(RayInfoContainer::const_iterator ri = rays.begin();
2588                ri != rays.end();
2589                ri++)
2590          if ((*ri).mRay->IsActive()) {
2591                Intersectable *object = (*ri).GetObject();
2592                if (object) {
2593                  if (!object->Mailed()) {
2594                        object->Mail();
2595                        float p = object->mCounter/(float)samples;
2596                        entropy -= p*log(p);
2597                  }
2598                }
2599          }
2600        entropy = entropy/log((float)samples);
2601  }
2602  else
2603        entropy = 1.0f;
2604 
2605  return entropy;
2606}
2607
2608float
2609RssTreeLeaf::ComputeRayLengthEntropy() const
2610{
2611  // get sum of all ray lengths
2612  // consider only passing rays or originating rays
2613  float sum = 0.0f;
2614  int samples = 0;
2615  int i=0;
2616  for(RayInfoContainer::const_iterator ri = rays.begin();
2617      ri != rays.end();
2618      ri++)
2619        if ((*ri).mRay->IsActive()) {
2620//        float s;
2621//        if (i == 0)
2622//              s = 200;
2623//        else
2624//              s = 100;
2625//        i++;
2626         
2627          sum += (*ri).mRay->GetSize();
2628          samples++;
2629        }
2630 
2631 
2632  float entropy = 0.0f;
2633 
2634  if (samples > 1) {
2635        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//              float p = s/sum;
2647                float p = (*ri).mRay->GetSize()/sum;
2648                entropy -= p*log(p);
2649          }
2650        entropy = entropy/log((float)samples);
2651  } else
2652        entropy = 1.0f;
2653 
2654  return entropy;
2655}
2656 
2657
2658
2659float
2660RssTreeLeaf::ComputeEntropyImportance() const
2661{
2662 
2663  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
2664  return 1.0f - ComputePvsEntropy();
2665}
2666
2667float
2668RssTreeLeaf::ComputeRayContributionImportance() const
2669{
2670  //  mPvsEntropy = ComputePvsEntropy();
2671  //  mRayLengthEntropy = ComputeRayLengthEntropy();
2672 
2673  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
2674  return 1.0f - ComputePvsEntropy();
2675}
2676
2677
2678AxisAlignedBox3
2679RssTree::GetShrankedBBox(const RssTreeNode *node) const
2680{
2681  if (node->parent == NULL)
2682        return bbox;
2683 
2684  if (!node->IsLeaf())
2685        return ((RssTreeInterior *)node)->bbox;
2686
2687  // evaluate bounding box from the ray origins
2688  AxisAlignedBox3 box;
2689  box.Initialize();
2690  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2691  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2692          ri != leaf->rays.end();
2693          ri++)
2694        if ((*ri).mRay->IsActive()) {
2695          box.Include((*ri).GetOrigin());
2696        }
2697  return box;
2698}
2699
2700
2701int
2702RssTree::CollectRays(VssRayContainer &rays,
2703                                         const int number)
2704{
2705  VssRayContainer allRays;
2706  CollectRays(allRays);
2707
2708 
2709  int desired = min(number, (int)allRays.size());
2710  float prob = desired/(float)allRays.size();
2711  while (rays.size() < desired) {
2712        VssRayContainer::const_iterator it = allRays.begin();
2713        for (; it != allRays.end() && rays.size() < desired; it++) {
2714          if (Random(1.0f) < prob)
2715                rays.push_back(*it);
2716        }
2717  }
2718  return (int)rays.size();
2719}
2720
2721
2722int
2723RssTree::CollectRays(VssRayContainer &rays
2724                                         )
2725{
2726  VssRay::NewMail();
2727
2728  stack<RssTreeNode *> tstack;
2729  PushRoots(tstack);
2730 
2731  while (!tstack.empty()) {
2732    RssTreeNode *node = tstack.top();
2733    tstack.pop();
2734       
2735    if (node->IsLeaf()) {
2736          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2737          // update pvs size
2738          RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
2739          for (;it != leaf->rays.end(); ++it)
2740                if (!(*it).mRay->Mailed()) {
2741                  (*it).mRay->Mail();
2742                  rays.push_back((*it).mRay);
2743                }
2744        } else {
2745          RssTreeInterior *in = (RssTreeInterior *)node;
2746          // both nodes for directional splits
2747          tstack.push(in->front);
2748          tstack.push(in->back);
2749        }
2750  }
2751 
2752  return (int)rays.size();
2753}
2754
2755
2756RssTreeNode *
2757RssTree::GetRoot(Intersectable *object) const
2758{
2759  if (mPerObjectTree && object) {
2760        int id = object->GetId();
2761        if (id >= mRoots.size())
2762          id = mRoots.size()-1; // $$ last tree is used by all unsigned objects
2763        return mRoots[id];
2764  } else
2765        return mRoots[0];
2766}
2767
2768void
2769RssTree::PushRoots(priority_queue<TraversalData> &st) const
2770{
2771  for (int i=0; i < mRoots.size(); i++)
2772        st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
2773}
2774
2775void
2776RssTree::PushRoots(stack<RssTreeNode *> &st)  const
2777{
2778  for (int i=0; i < mRoots.size(); i++)
2779        st.push(mRoots[i]);
2780}
2781
2782void
2783RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info)  const
2784{
2785  for (int i=0; i < mRoots.size(); i++)
2786        st.push(RayTraversalData(mRoots[i], info));
2787}
Note: See TracBrowser for help on using the repository browser.