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

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

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