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

Revision 492, 63.0 KB checked in by bittner, 19 years ago (diff)

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