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

Revision 500, 63.0 KB checked in by mattausch, 19 years ago (diff)

added new node primitive ro vspkd tree
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  return true; // corr. matt
1925}
1926
1927bool
1928RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
1929                                                                const RssTreeNode::RayInfo &r1,
1930                                                                const RssTreeNode::RayInfo &r2,
1931                                                                const RssTreeNode::RayInfo &r3)
1932{
1933 
1934
1935  return false;
1936}
1937
1938int
1939RssTree::RemoveInteriorRays(
1940                                                        RssTreeLeaf *leaf
1941                                                        )
1942{
1943#if 1
1944  // first collect all objects refered in this leaf
1945  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> > rayMap;
1946 
1947  RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
1948  for (; it != leaf->rays.end(); ++it) {
1949        Intersectable *object = (*it).GetObject();
1950       
1951        rayMap[object].push_back(&(*it));
1952//      vector<RayInfo *> *data = rayMap.Find(object);
1953//      if (data) {
1954//        data->push_back(&(*it));
1955//      } else {
1956//        //      rayMap[object] = vector<RayInfo *>;
1957//        rayMap[object].push_back(&(*it));
1958//      }
1959  }
1960
1961  // now go through all objects
1962  map<Intersectable *, vector<RssTreeLeaf::RayInfo *> >::iterator mi;
1963
1964  // infos of mailed rays are scheduled for removal
1965  VssRay::NewMail();
1966  for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
1967        vector<RssTreeLeaf::RayInfo *> &rays = (*mi).second;
1968        vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
1969        int rayIndex = 0;
1970       
1971        for (; ri != rays.end(); ++ri, ++rayIndex) {
1972          RssTreeNode::RayInfo *ray = *ri;
1973          int tries = rays.size();
1974          for (int i = 0; i < tries; i++) {
1975                int r1, r2, r3;
1976                GetRandomTripple(rays,
1977                                                 rayIndex,
1978                                                 r1,
1979                                                 r2,
1980                                                 r3);
1981                if (IsRayConvexCombination(*ray,
1982                                                                   *rays[r1],
1983                                                                   *rays[r2],
1984                                                                   *rays[r3])) {
1985                  ray->mRay->Mail();
1986                }
1987          }
1988        }
1989  }
1990 
1991
1992#endif 
1993        return 0;
1994}
1995
1996void
1997RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
1998                                                  const int numberOfRays,
1999                                                  SimpleRayContainer &rays)
2000{
2001
2002  if (numberOfRays == 0)
2003        return;
2004 
2005  int nrays = leaf->rays.size();
2006
2007
2008  AxisAlignedBox3 box = GetBBox(leaf);
2009  AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2010
2011  AxisAlignedBox3 sBox(box);
2012  AxisAlignedBox3 sDirBox(dirBox);
2013  const float smoothRange = 0.5f;
2014  sBox.Scale(1.0f + smoothRange);
2015  sDirBox.Scale(1.0f + smoothRange);
2016  int smoothRays = (int)numberOfRays*0.0f;
2017 
2018#if AVG_LEN
2019  float avgLength = 0.0f;
2020
2021  if (nrays) {
2022        const int numAvg = 5;
2023       
2024        int step = (nrays-1)/10;
2025        if (step<1)
2026          step = 1;
2027       
2028        int summed = 0;
2029        float sumLength = 0.0f;
2030        for (int i=0; i < nrays; i+=step) {
2031          sumLength += leaf->rays[i].mRay->GetSize();
2032          summed++;
2033        }
2034 
2035        avgLength = sumLength/summed;
2036  }
2037#endif
2038 
2039  Exporter *exporter = NULL;
2040  VssRayContainer selectedRays;
2041  static int counter = 0;
2042  if (counter < 0) {
2043        char filename[128];
2044        sprintf(filename, "selected-rays%04d.x3d", counter);
2045        exporter = Exporter::GetExporter(filename);
2046        //      exporter->SetWireframe();
2047        //      exporter->ExportKdTree(*mKdTree);
2048        exporter->SetWireframe();
2049        exporter->ExportScene(preprocessor->mSceneGraph->mRoot);
2050        exporter->SetWireframe();
2051        exporter->ExportBox(box);
2052        exporter->SetFilled();
2053        counter++;
2054  }
2055 
2056  int numberOfTries = numberOfRays*4;
2057  int generated = 0;
2058
2059  Intersectable::NewMail();
2060  vector<Intersectable *> pvs(0);
2061  pvs.reserve(leaf->GetPvsSize());
2062  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2063          ri != leaf->rays.end();
2064          ri++) {
2065        Intersectable *object = (*ri).GetObject();
2066        if (object && !object->Mailed()) {
2067          object->Mail();
2068          pvs.push_back(object);
2069          if (exporter)
2070                exporter->ExportIntersectable(object);
2071        }
2072  }
2073 
2074  for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2075        bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) < 0.7f));
2076       
2077        Vector3 origin, direction;
2078        // generate half of convex combination and half of random rays
2079        if (useExtendedConvexCombination) {
2080          // pickup 2 random rays
2081          int r1, r2;
2082          PickEdgeRays(leaf, r1, r2);
2083
2084         
2085          Vector3 o1 = leaf->rays[r1].GetOrigin();
2086          Vector3 o2 = leaf->rays[r2].GetOrigin();
2087         
2088          const float overlap = 0.0f;
2089          float w1, w2;
2090          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2091          origin = w1*o1 + w2*o2;
2092          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2093          direction =
2094                w1*leaf->rays[r1].GetDir() +
2095                w2*leaf->rays[r2].GetDir();
2096        } else {
2097          Vector3 dirVector;
2098
2099          Vector3 pVector;
2100          Vector3 dVector;
2101
2102          bool useHalton = true;
2103         
2104          if (useHalton) {
2105                // generate a random 5D vector
2106                pVector = Vector3(leaf->halton.GetNumber(1),
2107                                                  leaf->halton.GetNumber(2),
2108                                                  leaf->halton.GetNumber(3));
2109               
2110                dVector = Vector3(leaf->halton.GetNumber(4),
2111                                                  leaf->halton.GetNumber(5),
2112                                                  0.0f);
2113                leaf->halton.GenerateNext();
2114          } else {
2115                pVector = Vector3(RandomValue(0.0f, 1.0f),
2116                                                  RandomValue(0.0f, 1.0f),
2117                                                  RandomValue(0.0f, 1.0f));
2118               
2119                dVector = Vector3(RandomValue(0.0f, 1.0f),
2120                                                  RandomValue(0.0f, 1.0f),
2121                                                  0.0f);
2122          }
2123         
2124          if (i < smoothRays) {
2125                origin = sBox.GetPoint(pVector);
2126                dirVector = sDirBox.GetPoint(dVector);
2127          } else {
2128                origin = box.GetPoint(pVector);
2129                dirVector = dirBox.GetPoint(dVector);
2130          }
2131          direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
2132        }
2133
2134        // shift the origin a little bit
2135        direction.Normalize();
2136       
2137        //      float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2138        float dist = 0.0f;
2139        float minT, maxT;
2140        // compute interection of the ray with the box
2141
2142        if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2143          dist = maxT;
2144 
2145       
2146        origin += direction*dist;
2147       
2148        bool intersects = false;
2149       
2150
2151        if (i > numberOfRays/2) {
2152          if (exporter) {
2153                VssRay *ray = new VssRay(origin, origin + 100*direction, NULL, NULL);
2154                selectedRays.push_back(ray);
2155          }
2156
2157          // check whether the ray does not intersect already visible objects
2158          Ray traversalRay;
2159          traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2160          for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2161                  oi != pvs.end();
2162                  oi++) {
2163                Intersectable *object = *oi;
2164               
2165                if ( object->CastRay(traversalRay) ) {
2166                  intersects = true;
2167                  break;
2168                }
2169          }
2170        }
2171       
2172       
2173        if (!intersects) {
2174          //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2175          rays.push_back(SimpleRay(origin, direction));
2176          generated++;
2177        }
2178  }
2179
2180  //  cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2181 
2182  if (exporter) {
2183        exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2184        delete exporter;
2185        CLEAR_CONTAINER(selectedRays);
2186  }
2187}
2188
2189int
2190RssTree::PruneRays(RssTreeLeaf *leaf,
2191                                   const float contributionThreshold)
2192{
2193  int i;
2194  int j;
2195 
2196  for (j=0, i=0; i < leaf->rays.size(); i++) {
2197       
2198        if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
2199          // copy a valid sample
2200          if (i!=j)
2201                leaf->rays[j] = leaf->rays[i];
2202          j++;
2203        } else {
2204          // delete the ray
2205          leaf->rays[i].mRay->Unref();
2206          if (leaf->rays[i].mRay->RefCount() != 0) {
2207                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2208                exit(1);
2209          }
2210          delete leaf->rays[i].mRay;
2211        }
2212  }
2213
2214
2215  leaf->rays.resize(j);
2216  int removed = (i-j);
2217  stat.rayRefs -= removed;
2218  return removed;
2219}
2220
2221int
2222RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2223                                                 const float ratio)
2224{
2225  int i;
2226  int j;
2227 
2228  for (j=0, i=0; i < leaf->rays.size(); i++) {
2229       
2230        if (Random(1.0f) < ratio) {
2231          // copy a valid sample
2232          if (i!=j)
2233                leaf->rays[j] = leaf->rays[i];
2234          j++;
2235        } else {
2236          // delete the ray
2237          leaf->rays[i].mRay->Unref();
2238          if (leaf->rays[i].mRay->RefCount() != 0) {
2239                cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2240                exit(1);
2241          }
2242          delete leaf->rays[i].mRay;
2243        }
2244  }
2245
2246
2247  leaf->rays.resize(j);
2248  int removed = (i-j);
2249  stat.rayRefs -= removed;
2250  return removed;
2251}
2252
2253int
2254RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2255                                                           const float ratio)
2256{
2257  int i;
2258 
2259  if (leaf->rays.size() == 0)
2260        return 0;
2261 
2262  sort(leaf->rays.begin(),
2263           leaf->rays.end(),
2264           RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2265
2266  int desired = ratio*leaf->rays.size();
2267  int removed = leaf->rays.size() - desired;
2268 
2269  for (i=desired; i < leaf->rays.size(); i++) {
2270        // delete the ray
2271        leaf->rays[i].mRay->Unref();
2272        if (leaf->rays[i].mRay->RefCount() != 0) {
2273          cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2274          exit(1);
2275        }
2276        delete leaf->rays[i].mRay;
2277  }
2278 
2279  leaf->rays.resize(desired);
2280  stat.rayRefs -= removed;
2281  return removed;
2282}
2283
2284
2285int
2286RssTree::PruneRays(
2287                                   const int desired
2288                                   )
2289{
2290  bool globalPrunning = false;
2291
2292  stack<RssTreeNode *> tstack;
2293  int prunned = 0;
2294
2295  Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl;
2296
2297  if (globalPrunning) {
2298        VssRayContainer allRays;
2299        allRays.reserve(stat.rayRefs);
2300        CollectRays(allRays);
2301        sort(allRays.begin(),
2302                 allRays.end(),
2303                 GreaterWeightedPvsContribution);
2304       
2305        if ( desired >= allRays.size() )
2306          return 0;
2307       
2308        float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
2309       
2310        tstack.push(root);
2311       
2312        while (!tstack.empty()) {
2313          RssTreeNode *node = tstack.top();
2314          tstack.pop();
2315         
2316          if (node->IsLeaf()) {
2317                RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2318                prunned += PruneRays(leaf, contributionThreshold);
2319          } else {
2320                RssTreeInterior *in = (RssTreeInterior *)node;
2321                // both nodes for directional splits
2322                tstack.push(in->front);
2323                tstack.push(in->back);
2324          }
2325        }
2326  } else {
2327        // prune random rays from each leaf so that the ratio's remain the same
2328        tstack.push(root);
2329        float ratio = desired/(float)stat.rayRefs;
2330       
2331        while (!tstack.empty()) {
2332          RssTreeNode *node = tstack.top();
2333          tstack.pop();
2334         
2335          if (node->IsLeaf()) {
2336                RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2337                //              prunned += PruneRaysRandom(leaf, ratio);
2338                prunned += PruneRaysContribution(leaf, ratio);
2339          } else {
2340                RssTreeInterior *in = (RssTreeInterior *)node;
2341                // both nodes for directional splits
2342                tstack.push(in->front);
2343                tstack.push(in->back);
2344          }
2345        }
2346
2347       
2348
2349
2350  }
2351 
2352 
2353 
2354  Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl;
2355 
2356  return prunned;
2357}
2358
2359int
2360RssTree::GenerateRays(const int numberOfRays,
2361                                          const int numberOfLeaves,
2362                                          SimpleRayContainer &rays)
2363{
2364
2365
2366  vector<RssTreeLeaf *> leaves;
2367       
2368  CollectLeaves(leaves);
2369       
2370  sort(leaves.begin(),
2371           leaves.end(),
2372           GreaterContribution);
2373                         
2374
2375  float sumContrib = 0.0;
2376  int i;
2377  int k = 0;
2378  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2379        if (ValidLeaf(leaves[i])) {
2380          float c = leaves[i]->GetImportance();
2381          sumContrib += c;
2382          //            cout<<"ray contrib "<<i<<" : "<<c<<endl;
2383          k++;
2384        }
2385                               
2386  float avgContrib = sumContrib/numberOfLeaves;
2387 
2388  // always generate at leat n ray per leaf
2389  int fixedPerLeaf = 1;
2390  int fixed = 1*leaves.size();
2391  int iGenerated = numberOfRays;
2392  float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2393
2394  k = 0;
2395  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2396        if (ValidLeaf(leaves[i])) {
2397          k++;
2398          RssTreeLeaf *leaf = leaves[i];
2399          float c = leaf->GetImportance();
2400
2401          int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2402          GenerateLeafRays(leaf, num, rays);
2403        }
2404
2405  return rays.size();
2406}
2407
2408
2409float
2410RssTree::GetAvgPvsSize()
2411{
2412  stack<RssTreeNode *> tstack;
2413  tstack.push(root);
2414
2415  int sumPvs = 0;
2416  int leaves = 0;
2417  while (!tstack.empty()) {
2418    RssTreeNode *node = tstack.top();
2419    tstack.pop();
2420               
2421    if (node->IsLeaf()) {
2422          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2423          // update pvs size
2424          UpdatePvsSize(leaf);
2425          sumPvs += leaf->GetPvsSize();
2426          leaves++;
2427        } else {
2428          RssTreeInterior *in = (RssTreeInterior *)node;
2429          // both nodes for directional splits
2430          tstack.push(in->front);
2431          tstack.push(in->back);
2432        }
2433  }
2434
2435
2436  return sumPvs/(float)leaves;
2437}
2438
2439float weightAbsContributions = 0.0f;
2440// if small very high importance of the last sample
2441// if 1.0f then weighs = 1 1/2 1/3 1/4
2442float passSampleWeightDecay = 0.5f;
2443
2444void
2445RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2446                                                        float &weight,
2447                                                        float &contribution)
2448{
2449  VssRay *ray = info.mRay;
2450  weight = 1.0f/(mCurrentPass - ray->mPass + passSampleWeightDecay);
2451  contribution =
2452        weightAbsContributions*ray->mPvsContribution +
2453        (1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2454  // store the computed value
2455  info.mRay->mWeightedPvsContribution = weight*contribution;
2456}
2457
2458float
2459RssTree::GetSampleWeight(const int pass)
2460{
2461  int passDiff = mCurrentPass - pass;
2462  float weight;
2463  if (1)
2464        weight = 1.0f/(passDiff + passSampleWeightDecay);
2465  else
2466        switch (passDiff) {
2467        case 0:
2468          weight = 1.0f;
2469          break;
2470        default:
2471          weight = 0.01f;
2472          break;
2473          //    case 1:
2474//        weight = 0.5f;
2475//        break;
2476//      case 2:
2477//        weight = 0.25f;
2478//        break;
2479//      case 3:
2480//        weight = 0.12f;
2481//        break;
2482//      case 4:
2483//        weight = 0.06f;
2484//        break;
2485//      default:
2486//        weight = 0.03f;
2487//        break;
2488        }
2489  return weight;
2490}
2491
2492void
2493RssTree::ComputeImportance(RssTreeLeaf *leaf)
2494{
2495  if (0)
2496        leaf->mImportance = leaf->GetAvgRayContribution();
2497  else {
2498        RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2499          it_end = leaf->rays.end();
2500       
2501        float sumContributions = 0.0f;
2502        float sumRelContributions = 0.0f;
2503        float sumWeights = 0.0f;
2504        for (; it != it_end; ++it) {
2505          VssRay *ray = (*it).mRay;
2506          float weight = GetSampleWeight(ray->mPass);
2507         
2508          sumContributions += weight*ray->mPvsContribution;
2509          sumRelContributions += weight*ray->mRelativePvsContribution;
2510         
2511          //      sumWeights += weight;
2512          sumWeights += 1.0f;
2513        }
2514        // $$
2515        // sumWeights = leaf->mTotalRays;
2516         
2517        if (sumWeights != 0.0f)
2518          leaf->mImportance =
2519                (weightAbsContributions*sumContributions +
2520                 (1.0f - weightAbsContributions)*sumRelContributions)/sumWeights;
2521        else
2522          leaf->mImportance = 0.0f;
2523       
2524  }
2525 
2526  // return GetAvgRayContribution()*mImportance;
2527  //return GetAvgRayContribution();
2528}
2529
2530
2531float
2532RssTreeLeaf::GetImportance()  const
2533{
2534  return mImportance;
2535}
2536
2537
2538float
2539RssTreeLeaf::ComputePvsEntropy() const
2540{
2541  int samples = 0;
2542  Intersectable::NewMail();
2543  // set all object as belonging to the fron pvs
2544  for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2545          ri != rays.end();
2546          ri++)
2547        if ((*ri).mRay->IsActive()) {
2548          Intersectable *object = (*ri).GetObject();
2549          if (object) {
2550                if (!object->Mailed()) {
2551                  object->Mail();
2552                  object->mCounter = 1;
2553                } else
2554                  object->mCounter++;
2555                samples++;
2556          }
2557        }
2558 
2559  float entropy = 0.0f;
2560 
2561  if (samples > 1) {
2562        Intersectable::NewMail();
2563        for(RayInfoContainer::const_iterator ri = rays.begin();
2564                ri != rays.end();
2565                ri++)
2566          if ((*ri).mRay->IsActive()) {
2567                Intersectable *object = (*ri).GetObject();
2568                if (object) {
2569                  if (!object->Mailed()) {
2570                        object->Mail();
2571                        float p = object->mCounter/(float)samples;
2572                        entropy -= p*log(p);
2573                  }
2574                }
2575          }
2576        entropy = entropy/log((float)samples);
2577  }
2578  else
2579        entropy = 1.0f;
2580 
2581  return entropy;
2582}
2583
2584float
2585RssTreeLeaf::ComputeRayLengthEntropy() const
2586{
2587  // get sum of all ray lengths
2588  // consider only passing rays or originating rays
2589  float sum = 0.0f;
2590  int samples = 0;
2591  int i=0;
2592  for(RayInfoContainer::const_iterator ri = rays.begin();
2593      ri != rays.end();
2594      ri++)
2595        if ((*ri).mRay->IsActive()) {
2596//        float s;
2597//        if (i == 0)
2598//              s = 200;
2599//        else
2600//              s = 100;
2601//        i++;
2602         
2603          sum += (*ri).mRay->GetSize();
2604          samples++;
2605        }
2606 
2607 
2608  float entropy = 0.0f;
2609 
2610  if (samples > 1) {
2611        i = 0;
2612        for(RayInfoContainer::const_iterator ri = rays.begin();
2613                ri != rays.end();
2614                ri++)
2615          if ((*ri).mRay->IsActive()) {
2616//              float s;
2617//              if (i==0)
2618//                s = 200;
2619//              else
2620//                s = 100;
2621//              i++;
2622//              float p = s/sum;
2623                float p = (*ri).mRay->GetSize()/sum;
2624                entropy -= p*log(p);
2625          }
2626        entropy = entropy/log((float)samples);
2627  } else
2628        entropy = 1.0f;
2629 
2630  return entropy;
2631}
2632 
2633
2634
2635float
2636RssTreeLeaf::ComputeEntropyImportance() const
2637{
2638 
2639  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
2640  return 1.0f - ComputePvsEntropy();
2641}
2642
2643float
2644RssTreeLeaf::ComputeRayContributionImportance() const
2645{
2646  //  mPvsEntropy = ComputePvsEntropy();
2647  //  mRayLengthEntropy = ComputeRayLengthEntropy();
2648 
2649  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
2650  return 1.0f - ComputePvsEntropy();
2651}
2652
2653
2654AxisAlignedBox3
2655RssTree::GetShrankedBBox(const RssTreeNode *node) const
2656{
2657  if (node->parent == NULL)
2658        return bbox;
2659 
2660  if (!node->IsLeaf())
2661        return ((RssTreeInterior *)node)->bbox;
2662
2663  // evaluate bounding box from the ray origins
2664  AxisAlignedBox3 box;
2665  box.Initialize();
2666  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2667  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2668          ri != leaf->rays.end();
2669          ri++)
2670        if ((*ri).mRay->IsActive()) {
2671          box.Include((*ri).GetOrigin());
2672        }
2673  return box;
2674}
2675
2676
2677int
2678RssTree::CollectRays(VssRayContainer &rays,
2679                                         const int number)
2680{
2681  VssRayContainer allRays;
2682  CollectRays(allRays);
2683
2684 
2685  int desired = min(number, (int)allRays.size());
2686  float prob = desired/(float)allRays.size();
2687  while (rays.size() < desired) {
2688        VssRayContainer::const_iterator it = allRays.begin();
2689        for (; it != allRays.end() && rays.size() < desired; it++) {
2690          if (Random(1.0f) < prob)
2691                rays.push_back(*it);
2692        }
2693  }
2694  return (int)rays.size();
2695}
2696
2697
2698int
2699RssTree::CollectRays(VssRayContainer &rays
2700                                         )
2701{
2702  VssRay::NewMail();
2703
2704  stack<RssTreeNode *> tstack;
2705  tstack.push(root);
2706 
2707  while (!tstack.empty()) {
2708    RssTreeNode *node = tstack.top();
2709    tstack.pop();
2710       
2711    if (node->IsLeaf()) {
2712          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2713          // update pvs size
2714          RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
2715          for (;it != leaf->rays.end(); ++it)
2716                if (!(*it).mRay->Mailed()) {
2717                  (*it).mRay->Mail();
2718                  rays.push_back((*it).mRay);
2719                }
2720        } else {
2721          RssTreeInterior *in = (RssTreeInterior *)node;
2722          // both nodes for directional splits
2723          tstack.push(in->front);
2724          tstack.push(in->back);
2725        }
2726  }
2727 
2728  return (int)rays.size();
2729}
2730
Note: See TracBrowser for help on using the repository browser.