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

Revision 466, 48.1 KB checked in by bittner, 19 years ago (diff)

changed the viewcellsmanager interface to use vssrays - some functionality like the bsp merging is now restricted
Line 
1
2// ================================================================
3// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4// ****************************************************************
5/**
6   The KD tree based LSDS
7*/
8// Initial coding by
9/**
10   @author Jiri Bittner
11*/
12
13// Standard headers
14#include <stack>
15#include <queue>
16#include <algorithm>
17#include <fstream>
18#include <string>
19
20#include "RssTree.h"
21
22#include "Environment.h"
23#include "VssRay.h"
24#include "Intersectable.h"
25#include "Ray.h"
26#include "Containers.h"
27
28#define DEBUG_SPLIT_COST 0
29#define DEBUG_SPLITS 0
30
31// Static variables
32int
33RssTreeLeaf::mailID = 0;
34
35inline void
36AddObject2Pvs(Intersectable *object,
37                          const int side,
38                          int &pvsBack,
39                          int &pvsFront)
40{
41 
42  if (!object)
43        return;
44 
45  if (side <= 0) {
46        if (!object->Mailed() && !object->Mailed(2)) {
47          pvsBack++;
48          if (object->Mailed(1))
49                object->Mail(2);
50          else
51                object->Mail();
52        }
53  }
54 
55  if (side >= 0) {
56        if (!object->Mailed(1) && !object->Mailed(2)) {
57          pvsFront++;
58          if (object->Mailed())
59                object->Mail(2);
60          else
61                object->Mail(1);
62        }
63  }
64}
65
66// Constructor
67RssTree::RssTree()
68{
69  environment->GetIntValue("RssTree.maxDepth", termMaxDepth);
70  environment->GetIntValue("RssTree.minPvs", termMinPvs);
71  environment->GetIntValue("RssTree.minRays", termMinRays);
72  environment->GetFloatValue("RssTree.maxRayContribution", termMaxRayContribution);
73  environment->GetFloatValue("RssTree.maxCostRatio", termMaxCostRatio);
74
75  environment->GetFloatValue("RssTree.minSize", termMinSize);
76  termMinSize = sqr(termMinSize);
77       
78  environment->GetFloatValue("RssTree.refDirBoxMaxSize", refDirBoxMaxSize);
79  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
80 
81  environment->GetFloatValue("RssTree.epsilon", epsilon);
82  environment->GetFloatValue("RssTree.ct_div_ci", ct_div_ci);
83       
84  environment->GetFloatValue("RssTree.maxTotalMemory", maxTotalMemory);
85  environment->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
86 
87 
88
89
90  float refDirAngle;
91  environment->GetFloatValue("RssTree.refDirAngle", refDirAngle);
92 
93  environment->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
94  //= 1000;
95  environment->GetIntValue("RssTree.minCollapseDepth", minCollapseDepth);
96  //  int minCollapseDepth = 4;
97
98  //  pRefDirThresh = cos(0.5*M_PI - M_PI*refDirAngle/180.0);
99  //  cosRefDir = cos(M_PI*refDirAngle/180.0);
100  //  sinRefDir = sin(M_PI*refDirAngle/180.0);
101 
102 
103  // split type
104  char sname[128];
105  environment->GetStringValue("RssTree.splitType", sname);
106  string name(sname);
107       
108  if (name.compare("regular") == 0)
109    splitType = ESplitRegular;
110  else
111    if (name.compare("heuristic") == 0)
112      splitType = ESplitHeuristic;
113        else
114          if (name.compare("hybrid") == 0)
115                splitType = ESplitHybrid;
116          else {
117                cerr<<"Invalid RssTree split type "<<name<<endl;
118                exit(1);
119          }
120       
121  environment->GetBoolValue("RssTree.randomize", randomize);
122  environment->GetBoolValue("RssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
123
124  environment->GetBoolValue("RssTree.interleaveDirSplits", mInterleaveDirSplits);
125  environment->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);
126
127  root = NULL;
128 
129  splitCandidates = new vector<SortableEntry>;
130}
131
132
133RssTree::~RssTree()
134{
135  if (root)
136    delete root;
137}
138
139
140
141
142void
143RssStatistics::Print(ostream &app) const
144{
145  app << "===== RssTree statistics ===============\n";
146
147  app << "#N_RAYS ( Number of rays )\n"
148      << rays <<endl;
149
150  app << "#N_INITPVS ( Initial PVS size )\n"
151      << initialPvsSize <<endl;
152 
153  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
154
155  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
156
157  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
158  for (int i=0; i<7; i++)
159    app << splits[i] <<" ";
160  app <<endl;
161
162  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
163    rayRefs << "\n";
164
165  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
166    rayRefs/(double)rays << "\n";
167
168  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
169    rayRefs/(double)Leaves() << "\n";
170
171  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
172    maxRayRefs << "\n";
173
174
175  //  app << setprecision(4);
176
177  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
178    maxDepthNodes*100/(double)Leaves()<<endl;
179
180  app << "#N_PMINPVSLEAVES  ( Percentage of leaves with minCost )\n"<<
181    minPvsNodes*100/(double)Leaves()<<endl;
182
183  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\n"<<
184    minRaysNodes*100/(double)Leaves()<<endl;
185       
186  app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
187    minSizeNodes*100/(double)Leaves()<<endl;
188
189  app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
190    maxRayContribNodes*100/(double)Leaves()<<endl;
191
192  app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
193    maxCostRatioNodes*100/(double)Leaves()<<endl;
194
195  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
196    addedRayRefs<<endl;
197
198  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
199    removedRayRefs<<endl;
200
201  //  app << setprecision(4);
202
203  app << "#N_CTIME  ( Construction time [s] )\n"
204      << Time() << " \n";
205
206  app << "===== END OF RssTree statistics ==========\n";
207
208}
209
210
211void
212RssTreeLeaf::UpdatePvsSize()
213{
214  if (!mValidPvs) {
215        Intersectable::NewMail();
216        int pvsSize = 0;
217        for(RssTreeNode::RayInfoContainer::iterator ri = rays.begin();
218                ri != rays.end();
219                ri++)
220          if ((*ri).mRay->IsActive()) {
221                Intersectable *object;
222#if BIDIRECTIONAL_RAY
223                object = (*ri).mRay->mOriginObject;
224                if (object && !object->Mailed()) {
225                  pvsSize++;
226                  object->Mail();
227                }
228#endif
229                object = (*ri).mRay->mTerminationObject;
230                if (object && !object->Mailed()) {
231                  pvsSize++;
232                  object->Mail();
233                }
234          }
235        mPvsSize = pvsSize;
236        mValidPvs = true;
237
238        ComputeEntropyImportance();
239  }
240}
241
242bool
243RssTree::ClipRay(
244                                 RssTreeNode::RayInfo &rayInfo,
245                                 const AxisAlignedBox3 &box
246                                 )
247{
248  float tmin, tmax;
249  static Ray ray;
250  ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
251  box.ComputeMinMaxT(ray, &tmin, &tmax);
252  if (tmin >= tmax)
253        return false;
254 
255  // now check if the ray origin lies inside the box
256  if ( tmax < rayInfo.mRay->GetSize() ) {
257        // this ray does not leave the box
258        rayInfo.mRay->SetupEndPoints(
259                                                                 ray.Extrap(tmax),
260                                                                 rayInfo.mRay->mTermination
261                                                                 );
262        return true;
263  }
264
265  return false;
266}
267
268
269void
270RssTree::Construct(
271                                   VssRayContainer &rays,
272                                   // forced bounding box is only used when computing from-box
273                                   // visibility
274                                   AxisAlignedBox3 *forcedBoundingBox
275                                   )
276{
277  stat.Start();
278 
279  maxMemory = maxStaticMemory;
280
281  if (root)
282    delete root;
283       
284  root = new RssTreeLeaf(NULL, rays.size());
285  // first construct a leaf that will get subdivide
286  RssTreeLeaf *leaf = (RssTreeLeaf *) root;
287
288  stat.nodes = 1;
289 
290  bbox.Initialize();
291  dirBBox.Initialize();
292
293
294  mForcedBoundingBox = forcedBoundingBox;
295  for(VssRayContainer::const_iterator ri = rays.begin();
296      ri != rays.end();
297      ri++) {
298               
299        RssTreeNode::RayInfo info(*ri);
300        if (forcedBoundingBox)
301          if (!ClipRay(info, *forcedBoundingBox))
302                continue;
303
304        leaf->AddRay(info);
305               
306    bbox.Include((*ri)->GetOrigin());
307    bbox.Include((*ri)->GetTermination());
308   
309               
310        dirBBox.Include(Vector3(
311                                                        (*ri)->GetDirParametrization(0),
312                                                        (*ri)->GetDirParametrization(1),
313                                                        0
314                                                        )
315                                        );
316  }
317       
318       
319  if ( forcedBoundingBox )
320        bbox = *forcedBoundingBox;
321       
322  cout<<"Bbox = "<<bbox<<endl;
323  cout<<"Dirr Bbox = "<<dirBBox<<endl;
324
325  stat.rays = leaf->rays.size();
326  leaf->UpdatePvsSize();
327  leaf->ComputeEntropyImportance();
328
329  stat.initialPvsSize = leaf->GetPvsSize();
330  // Subdivide();
331  root = Subdivide(TraversalData(leaf, bbox, 0));
332
333  if (splitCandidates) {
334    // force realease of this vector
335    delete splitCandidates;
336    splitCandidates = new vector<SortableEntry>;
337  }
338 
339  stat.Stop();
340
341  stat.Print(cout);
342  cout<<"#Total memory="<<GetMemUsage()<<endl;
343
344}
345
346int
347RssTree::UpdateSubdivision()
348{
349  priority_queue<TraversalData> tStack;
350  //  stack<TraversalData> tStack;
351 
352  tStack.push(TraversalData(root, bbox, 0));
353       
354  AxisAlignedBox3 backBox;
355  AxisAlignedBox3 frontBox;
356
357  maxMemory = maxTotalMemory;
358  int subdivided = 0;
359  int lastMem = 0;
360  while (!tStack.empty()) {
361               
362        float mem = GetMemUsage();
363               
364        if ( lastMem/10 != ((int)mem)/10) {
365          cout<<mem<<" MB"<<endl;
366        }
367        lastMem = (int)mem;
368               
369        if (  mem > maxMemory ) {
370      // count statistics on unprocessed leafs
371      while (!tStack.empty()) {
372                //                              EvaluateLeafStats(tStack.top());
373                tStack.pop();
374      }
375      break;
376    }
377   
378    TraversalData data = tStack.top();
379    tStack.pop();
380
381        if (data.node->IsLeaf()) {
382          RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
383                                                                                data.bbox,
384                                                                                backBox,
385                                                                                frontBox
386                                                                                );
387          if (!node->IsLeaf()) {
388                subdivided++;
389                RssTreeInterior *interior = (RssTreeInterior *) node;
390                // push the children on the stack
391                tStack.push(TraversalData(interior->back, backBox, data.depth+1));
392                tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
393          } else {
394                //      EvaluateLeafStats(data);
395          }
396        } else {
397          RssTreeInterior *interior = (RssTreeInterior *) data.node;
398          tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
399          tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
400        }
401  }
402  return subdivided;
403}
404
405
406RssTreeNode *
407RssTree::Subdivide(const TraversalData &tdata)
408{
409  RssTreeNode *result = NULL;
410
411  priority_queue<TraversalData> tStack;
412  //  stack<TraversalData> tStack;
413 
414  tStack.push(tdata);
415
416  AxisAlignedBox3 backBox;
417  AxisAlignedBox3 frontBox;
418
419 
420  int lastMem = 0;
421  while (!tStack.empty()) {
422
423        float mem = GetMemUsage();
424               
425        if ( lastMem/10 != ((int)mem)/10) {
426          cout<<mem<<" MB"<<endl;
427        }
428        lastMem = (int)mem;
429               
430        if (  mem > maxMemory ) {
431      // count statistics on unprocessed leafs
432      while (!tStack.empty()) {
433                EvaluateLeafStats(tStack.top());
434                tStack.pop();
435      }
436      break;
437    }
438   
439    TraversalData data = tStack.top();
440    tStack.pop();
441
442#if DEBUG_SPLITS
443        Debug<<"#Splitting node"<<endl;
444        data.node->Print(Debug);
445#endif
446        RssTreeNode *node = SubdivideNode((RssTreeLeaf *) data.node,
447                                                                          data.bbox,
448                                                                          backBox,
449                                                                          frontBox
450                                                                          );
451       
452
453        if (result == NULL)
454      result = node;
455   
456    if (!node->IsLeaf()) {
457                       
458      RssTreeInterior *interior = (RssTreeInterior *) node;
459      // push the children on the stack
460      tStack.push(TraversalData(interior->back, backBox, data.depth+1));
461      tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
462
463#if DEBUG_SPLITS
464          Debug<<"#New nodes"<<endl;
465          interior->back->Print(Debug);
466          interior->front->Print(Debug);
467          Debug<<"#####################################"<<endl;
468#endif
469         
470    } else {
471      EvaluateLeafStats(data);
472    }
473  }
474
475  return result;
476}
477
478
479// returns selected plane for subdivision
480int
481RssTree::SelectPlane(
482                                         RssTreeLeaf *leaf,
483                                         const AxisAlignedBox3 &box,
484                                         float &position,
485                                         int &raysBack,
486                                         int &raysFront,
487                                         int &pvsBack,
488                                         int &pvsFront
489                                         )
490{
491
492  int minDirDepth = 6;
493  int axis;
494  float costRatio;
495       
496  costRatio = BestCostRatio(leaf,
497                                                        axis,
498                                                        position,
499                                                        raysBack,
500                                                        raysFront,
501                                                        pvsBack,
502                                                        pvsFront
503                                                        );
504#if DEBUG_SPLIT_COST
505  cout<<axis<<" r="<<costRatio<<endl;
506#endif
507  if (costRatio > termMaxCostRatio) {
508        //              cout<<"Too big cost ratio "<<costRatio<<endl;
509        stat.maxCostRatioNodes++;
510        return -1;
511  }
512       
513#if 0   
514  cout<<
515        "pvs="<<leaf->mPvsSize<<
516        " rays="<<leaf->rays.size()<<
517        " rc="<<leaf->GetAvgRayContribution()<<
518        " axis="<<axis<<endl;
519#endif
520       
521  return axis;
522}
523
524
525float
526RssTree::GetCostRatio(
527                                          RssTreeLeaf *leaf,
528                                          const int axis,
529                                          const float position,
530                                          const int raysBack,
531                                          const int raysFront,
532                                          const int pvsBack,
533                                          const int pvsFront
534                                          )
535{
536  bool costImportance = true;
537
538  float ratio;
539  AxisAlignedBox3 box;
540  float minBox, maxBox;
541
542  if (axis < 3) {
543        box = GetBBox(leaf);
544        minBox = box.Min(axis);
545        maxBox = box.Max(axis);
546  }     else {
547        box = GetDirBBox(leaf);
548        minBox = box.Min(axis-3);
549        maxBox = box.Max(axis-3);
550  }
551       
552  float sizeBox = maxBox - minBox;
553
554  int pvsSize = leaf->GetPvsSize();
555
556  if (!costImportance) {
557        //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
558        float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
559        float newCost = ct_div_ci + sum/sizeBox;
560        float oldCost = pvsSize;
561        ratio = newCost/oldCost;
562  } else {
563        // importance based cost
564#if 0
565        float newContrib =
566          ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
567           (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
568               
569        //                      float newContrib =
570        //                              sqr(pvsBack/(raysBack + Limits::Small)) +
571        //                              sqr(pvsFront/(raysFront + Limits::Small));
572        float oldContrib = sqr(leaf->GetAvgRayContribution());
573        ratio = oldContrib/newContrib;
574#else
575#if 1
576        float newCost = raysBack*pvsBack  + raysFront*pvsFront;
577        float oldCost = leaf->rays.size()*pvsSize;
578        ratio = newCost/oldCost;
579#else
580#if 0
581        float newCost = (pvsBack  + pvsFront)*0.5f;
582        float oldCost = pvsSize;
583        ratio = newCost/oldCost;
584#else
585        float newCost = abs(raysBack  - raysFront);
586        float oldCost = leaf->rays.size();
587        ratio = newCost/oldCost;
588#endif
589#endif
590#endif
591  }
592       
593  return ratio;
594}
595                                                       
596
597float
598RssTree::EvalCostRatio(
599                                           RssTreeLeaf *leaf,
600                                           const int axis,
601                                           const float position,
602                                           int &raysBack,
603                                           int &raysFront,
604                                           int &pvsBack,
605                                           int &pvsFront
606                                           )
607{
608  raysBack = 0;
609  raysFront = 0;
610  pvsFront = 0;
611  pvsBack = 0;
612
613       
614  Intersectable::NewMail(3);
615       
616  if (axis <= RssTreeNode::SPLIT_Z) {
617    // this is the main ray classification loop!
618    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
619                ri != leaf->rays.end();
620                ri++)
621      if ((*ri).mRay->IsActive()) {
622                               
623                // determine the side of this ray with respect to the plane
624                int side = (*ri).ComputeRaySide(axis, position);
625                //                              (*ri).mRay->mSide = side;
626                               
627                if (side <= 0)
628                  raysBack++;
629                               
630                if (side >= 0)
631                  raysFront++;
632                               
633                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
634      }
635               
636  } else {
637               
638        // directional split
639    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
640                ri != leaf->rays.end();
641                ri++)
642      if ((*ri).mRay->IsActive()) {
643                               
644                // determine the side of this ray with respect to the plane
645                int side;
646                if ((*ri).mRay->GetDirParametrization(axis - 3) <= position)
647                  side = -1;
648                else
649                  side = 1;
650               
651                if (side <= 0)
652                  raysBack++;
653                               
654                if (side >= 0)
655                  raysFront++;
656
657                //                              (*ri).mRay->mSide = side;
658                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
659                               
660      }
661  }
662
663  float ratio = GetCostRatio(
664                                                         leaf,
665                                                         axis,
666                                                         position,
667                                                         raysBack,
668                                                         raysFront,
669                                                         pvsBack,
670                                                         pvsFront);
671
672  //    cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
673  //  float oldCost = leaf->rays.size();
674
675  //    cout<<"ratio="<<ratio<<endl;
676
677  return ratio;
678}
679
680float
681RssTree::BestCostRatio(
682                                           RssTreeLeaf *leaf,
683                                           int &axis,
684                                           float &position,
685                                           int &raysBack,
686                                           int &raysFront,
687                                           int &pvsBack,
688                                           int &pvsFront
689                                           )
690{
691  int nRaysBack[6], nRaysFront[6];
692  int nPvsBack[6], nPvsFront[6];
693  float nPosition[6];
694  float nCostRatio[6];
695  int bestAxis = -1;
696       
697  AxisAlignedBox3 sBox = GetBBox(leaf);
698  AxisAlignedBox3 dBox = GetDirBBox(leaf);
699  // int sAxis = box.Size().DrivingAxis();
700  int sAxis = sBox.Size().DrivingAxis();
701  int dAxis = dBox.Size().DrivingAxis() + 3;
702
703
704  float dirSplitBoxSize = 0.01f;
705  bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
706               
707 
708  for (axis = 0; axis < 5; axis++)
709        if (
710                (axis < 3 && (leaf->depth < mDirSplitDepth ||  mInterleaveDirSplits)) ||
711                (axis >= 3 && (leaf->depth >= mDirSplitDepth))
712                ) {
713          if (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis) {
714               
715                if (splitType == ESplitRegular) {
716                  if (axis < 3)
717                        nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
718                  else
719                        nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
720                 
721                  nCostRatio[axis] = EvalCostRatio(leaf,
722                                                                                   axis,
723                                                                                   nPosition[axis],
724                                                                                   nRaysBack[axis],
725                                                                                   nRaysFront[axis],
726                                                                                   nPvsBack[axis],
727                                                                                   nPvsFront[axis]
728                                                                                   );
729                } else
730                  if (splitType == ESplitHeuristic) {
731                        nCostRatio[axis] = EvalCostRatioHeuristic(
732                                                                                                          leaf,
733                                                                                                          axis,
734                                                                                                          nPosition[axis],
735                                                                                                          nRaysBack[axis],
736                                                                                                          nRaysFront[axis],
737                                                                                                          nPvsBack[axis],
738                                                                                                          nPvsFront[axis]);
739                  } else
740                        if (splitType == ESplitHybrid) {
741                          if (leaf->depth > 7)
742                                nCostRatio[axis] = EvalCostRatioHeuristic(
743                                                                                                                  leaf,
744                                                                                                                  axis,
745                                                                                                                  nPosition[axis],
746                                                                                                                  nRaysBack[axis],
747                                                                                                                  nRaysFront[axis],
748                                                                                                                  nPvsBack[axis],
749                                                                                                                  nPvsFront[axis]);
750                          else {
751                                if (axis < 3)
752                                  nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
753                                else
754                                  nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
755                               
756                                nCostRatio[axis] = EvalCostRatio(leaf,
757                                                                                                 axis,
758                                                                                                 nPosition[axis],
759                                                                                                 nRaysBack[axis],
760                                                                                                 nRaysFront[axis],
761                                                                                                 nPvsBack[axis],
762                                                                                                 nPvsFront[axis]
763                                                                                                 );
764                          }
765                        } else {
766                          cerr<<"RssTree: Unknown split heuristics\n";
767                          exit(1);
768                        }
769               
770               
771                if ( bestAxis == -1)
772                  bestAxis = axis;
773                else
774                  if ( nCostRatio[axis] < nCostRatio[bestAxis] )
775                        bestAxis = axis;
776          }
777        }
778 
779  axis = bestAxis;
780  position = nPosition[bestAxis];
781
782  raysBack = nRaysBack[bestAxis];
783  raysFront = nRaysFront[bestAxis];
784
785  pvsBack = nPvsBack[bestAxis];
786  pvsFront = nPvsFront[bestAxis];
787       
788  return nCostRatio[bestAxis];
789}
790
791       
792float
793RssTree::EvalCostRatioHeuristic(
794                                                                RssTreeLeaf *leaf,
795                                                                const int axis,
796                                                                float &bestPosition,
797                                                                int &raysBack,
798                                                                int &raysFront,
799                                                                int &pvsBack,
800                                                                int &pvsFront
801                                                                )
802{
803  AxisAlignedBox3 box;
804  float minBox, maxBox;
805       
806  if (axis < 3) {
807        box = GetBBox(leaf);
808        minBox = box.Min(axis);
809        maxBox = box.Max(axis);
810  } else {
811        box = GetDirBBox(leaf);
812        minBox = box.Min(axis-3);
813        maxBox = box.Max(axis-3);
814  }
815       
816  SortSplitCandidates(leaf, axis);
817 
818  // go through the lists, count the number of objects left and right
819  // and evaluate the following cost funcion:
820  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
821       
822  int rl=0, rr = leaf->rays.size();
823  int pl=0, pr = leaf->GetPvsSize();
824  float sizeBox = maxBox - minBox;
825       
826  float minBand = minBox + 0.1*(maxBox - minBox);
827  float maxBand = minBox + 0.9*(maxBox - minBox);
828       
829  float minRatio = 1e20;
830       
831  Intersectable::NewMail();
832  // set all object as belonging to the fron pvs
833  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
834          ri != leaf->rays.end();
835          ri++)
836        if ((*ri).mRay->IsActive()) {
837          Intersectable *object = (*ri).mRay->mTerminationObject;
838          if (object)
839                if (!object->Mailed()) {
840                  object->Mail();
841                  object->mCounter = 1;
842                } else
843                  object->mCounter++;
844        }
845       
846  Intersectable::NewMail();
847       
848  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
849          ci < splitCandidates->end();
850          ci++) {
851        VssRay *ray;
852        switch ((*ci).type) {
853        case SortableEntry::ERayMin: {
854          rr--;
855          rl++;
856          ray = (VssRay *) (*ci).data;
857          Intersectable *object = ray->mTerminationObject;
858          if (object) {
859                if (!object->Mailed()) {
860                  object->Mail();
861                  pl++;
862                }
863                if (--object->mCounter == 0)
864                  pr--;
865          }
866          break;
867        }
868        }
869
870        float position = (*ci).value;
871       
872        if (position > minBand && position < maxBand) {
873                       
874          float ratio = GetCostRatio(
875                                                                 leaf,
876                                                                 axis,
877                                                                 position,
878                                                                 rl,
879                                                                 rr,
880                                                                 pl,
881                                                                 pr);
882                       
883                       
884          //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
885          //      cout<<"cost= "<<sum<<endl;
886                       
887          if (ratio < minRatio) {
888                minRatio = ratio;
889                bestPosition = position;
890                               
891                raysBack = rl;
892                raysFront = rr;
893                               
894                pvsBack = pl;
895                pvsFront = pr;
896                               
897      }
898    }
899  }
900
901 
902  //  cout<<"===================="<<endl;
903  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
904  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
905  return minRatio;
906}
907
908void
909RssTree::SortSplitCandidates(
910                                                         RssTreeLeaf *node,
911                                                         const int axis
912                                                         )
913{
914 
915  splitCandidates->clear();
916 
917  int requestedSize = 2*(node->rays.size());
918  // creates a sorted split candidates array
919  if (splitCandidates->capacity() > 500000 &&
920      requestedSize < (int)(splitCandidates->capacity()/10) ) {
921   
922    delete splitCandidates;
923    splitCandidates = new vector<SortableEntry>;
924  }
925 
926  splitCandidates->reserve(requestedSize);
927
928  // insert all queries
929  for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
930      ri < node->rays.end();
931      ri++) {
932        if ((*ri).mRay->IsActive()) {
933          if (axis < 3) {
934                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
935                                                                                                 (*ri).ExtrapOrigin(axis),
936                                                                                                 (void *)(*ri).mRay)
937                                                                   );
938          } else {
939                float pos = (*ri).mRay->GetDirParametrization(axis-3);
940                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
941                                                                                                 pos,
942                                                                                                 (void *)(*ri).mRay)
943                                                                   );
944          }
945        }
946  }
947 
948  stable_sort(splitCandidates->begin(), splitCandidates->end());
949}
950
951
952void
953RssTree::EvaluateLeafStats(const TraversalData &data)
954{
955
956  // the node became a leaf -> evaluate stats for leafs
957  RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
958
959  if (data.depth >= termMaxDepth)
960    stat.maxDepthNodes++;
961 
962  //  if ( (int)(leaf->rays.size()) < termMinCost)
963  //    stat.minCostNodes++;
964  if ( leaf->GetPvsSize() < termMinPvs)
965        stat.minPvsNodes++;
966
967  if ( leaf->GetPvsSize() < termMinRays)
968        stat.minRaysNodes++;
969
970  if (leaf->GetAvgRayContribution() > termMaxRayContribution )
971        stat.maxRayContribNodes++;
972       
973  if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
974        stat.minSizeNodes++;
975  }
976
977  if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
978    stat.maxRayRefs = leaf->rays.size();
979
980}
981
982bool
983RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
984{
985  return ( (leaf->GetPvsSize() < termMinPvs) ||
986                   (leaf->rays.size() < termMinRays) ||
987                   //                    (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
988                   (leaf->depth >= termMaxDepth) ||
989                   (SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
990                   );
991}
992
993
994RssTreeNode *
995RssTree::SubdivideNode(
996                                           RssTreeLeaf *leaf,
997                                           const AxisAlignedBox3 &box,
998                                           AxisAlignedBox3 &backBBox,
999                                           AxisAlignedBox3 &frontBBox
1000                                           )
1001{
1002 
1003  if (TerminationCriteriaSatisfied(leaf)) {
1004#if 0
1005        if (leaf->depth >= termMaxDepth) {
1006          cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
1007          cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
1008        }
1009#endif
1010               
1011        return leaf;
1012  }
1013       
1014  float position;
1015       
1016  // first count ray sides
1017  int raysBack;
1018  int raysFront;
1019  int pvsBack;
1020  int pvsFront;
1021       
1022  // select subdivision axis
1023  int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
1024  Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
1025 
1026  if (axis == -1) {
1027    return leaf;
1028  }
1029 
1030  stat.nodes+=2;
1031  stat.splits[axis]++;
1032
1033  // add the new nodes to the tree
1034  RssTreeInterior *node = new RssTreeInterior(leaf->parent);
1035
1036  node->axis = axis;
1037  node->position = position;
1038  node->bbox = box;
1039  node->dirBBox = GetDirBBox(leaf);
1040 
1041  backBBox = box;
1042  frontBBox = box;
1043
1044  RssTreeLeaf *back = new RssTreeLeaf(node, raysBack);
1045  RssTreeLeaf *front = new RssTreeLeaf(node, raysFront);
1046
1047  // replace a link from node's parent
1048  if (  leaf->parent )
1049    leaf->parent->ReplaceChildLink(leaf, node);
1050  // and setup child links
1051  node->SetupChildLinks(back, front);
1052       
1053  if (axis <= RssTreeNode::SPLIT_Z) {
1054        backBBox.SetMax(axis, position);
1055    frontBBox.SetMin(axis, position);
1056               
1057        for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1058                ri != leaf->rays.end();
1059                ri++) {
1060      if ((*ri).mRay->IsActive()) {
1061                               
1062                // first unref ray from the former leaf
1063                (*ri).mRay->Unref();
1064
1065                // Debug << "computed t: " << (*ri).mRay->mT << endl;
1066                // determine the side of this ray with respect to the plane
1067                int side = node->ComputeRaySide(*ri);
1068               
1069                if (side == 0) {
1070                  if ((*ri).mRay->HasPosDir(axis)) {
1071                        back->AddRay(*ri);
1072                        front->AddRay(*ri);
1073                  } else {
1074                        back->AddRay(*ri);
1075                        front->AddRay(*ri);
1076                  }
1077                } else
1078                  if (side == 1)
1079                        front->AddRay(*ri);
1080                  else
1081                        back->AddRay(*ri);
1082      } else
1083                (*ri).mRay->Unref();
1084    }
1085  } else {
1086    // rays front/back
1087   
1088   
1089    for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1090                ri != leaf->rays.end();
1091                ri++) {
1092      if ((*ri).mRay->IsActive()) {
1093                // first unref ray from the former leaf
1094                (*ri).mRay->Unref();
1095
1096                int side;
1097                if ((*ri).mRay->GetDirParametrization(axis - 3) <= position)
1098                  side = -1;
1099                else
1100                  side = 1;
1101                               
1102                if (side == 1)
1103                  front->AddRay(*ri);
1104                else
1105                  back->AddRay(*ri);
1106                               
1107      } else
1108                (*ri).mRay->Unref();
1109    }
1110  }
1111
1112#if 0
1113  front->SetPvsSize(pvsFront);
1114  back->SetPvsSize(pvsBack);
1115  // compute entropy as well
1116  front->ComputeEntropyImportance();
1117  back->ComputeEntropyImportance();
1118#else
1119  front->UpdatePvsSize();
1120  back->UpdatePvsSize();
1121#endif
1122
1123 
1124  // update stats
1125  stat.rayRefs -= leaf->rays.size();
1126  stat.rayRefs += raysBack + raysFront;
1127
1128 
1129  delete leaf;
1130  return node;
1131}
1132
1133
1134
1135
1136
1137
1138int
1139RssTree::ReleaseMemory(const int time)
1140{
1141  stack<RssTreeNode *> tstack;
1142 
1143  // find a node in the tree which subtree will be collapsed
1144  int maxAccessTime = time - accessTimeThreshold;
1145  int released;
1146
1147  tstack.push(root);
1148
1149  while (!tstack.empty()) {
1150    RssTreeNode *node = tstack.top();
1151    tstack.pop();
1152   
1153 
1154    if (!node->IsLeaf()) {
1155      RssTreeInterior *in = (RssTreeInterior *)node;
1156      //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1157      if (in->depth >= minCollapseDepth &&
1158                  in->lastAccessTime <= maxAccessTime) {
1159                released = CollapseSubtree(node, time);
1160                break;
1161      }
1162     
1163      if (in->back->GetAccessTime() < 
1164                  in->front->GetAccessTime()) {
1165                tstack.push(in->front);
1166                tstack.push(in->back);
1167      } else {
1168                tstack.push(in->back);
1169                tstack.push(in->front);
1170      }
1171    }
1172  }
1173
1174  while (tstack.empty()) {
1175    // could find node to collaps...
1176    //    cout<<"Could not find a node to release "<<endl;
1177    break;
1178  }
1179 
1180  return released;
1181}
1182
1183
1184
1185
1186RssTreeNode *
1187RssTree::SubdivideLeaf(
1188                                           RssTreeLeaf *leaf
1189                                           )
1190{
1191  RssTreeNode *node = leaf;
1192       
1193  AxisAlignedBox3 leafBBox = GetBBox(leaf);
1194
1195  static int pass = 0;
1196  pass ++;
1197       
1198  // check if we should perform a dynamic subdivision of the leaf
1199  if (!TerminationCriteriaSatisfied(leaf)) {
1200   
1201        // memory check and realese...
1202    if (GetMemUsage() > maxTotalMemory) {
1203      ReleaseMemory( pass );
1204    }
1205   
1206    AxisAlignedBox3 backBBox, frontBBox;
1207
1208    // subdivide the node
1209    node =
1210      SubdivideNode(leaf,
1211                                        leafBBox,
1212                                        backBBox,
1213                                        frontBBox
1214                                        );
1215  }
1216       
1217  return node;
1218}
1219
1220
1221
1222void
1223RssTree::UpdateRays(VssRayContainer &remove,
1224                                        VssRayContainer &add
1225                                        )
1226{
1227  RssTreeLeaf::NewMail();
1228
1229  // schedule rays for removal
1230  for(VssRayContainer::const_iterator ri = remove.begin();
1231      ri != remove.end();
1232      ri++) {
1233    (*ri)->ScheduleForRemoval();
1234  }
1235
1236  int inactive=0;
1237
1238  for(VssRayContainer::const_iterator ri = remove.begin();
1239      ri != remove.end();
1240      ri++) {
1241    if ((*ri)->ScheduledForRemoval())
1242          //      RemoveRay(*ri, NULL, false);
1243          // !!! BUG - with true it does not work correctly - aggreated delete
1244      RemoveRay(*ri, NULL, true);
1245    else
1246      inactive++;
1247  }
1248
1249
1250  //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1251 
1252  for(VssRayContainer::const_iterator ri = add.begin();
1253      ri != add.end();
1254      ri++) {
1255        RssTreeNode::RayInfo info(*ri);
1256        if (mForcedBoundingBox==NULL || ClipRay(info, bbox))
1257          AddRay(info);
1258  }
1259}
1260
1261
1262void
1263RssTree::RemoveRay(VssRay *ray,
1264                                   vector<RssTreeLeaf *> *affectedLeaves,
1265                                   const bool removeAllScheduledRays
1266                                   )
1267{
1268       
1269  stack<RayTraversalData> tstack;
1270
1271  tstack.push(RayTraversalData(root, RssTreeNode::RayInfo(ray)));
1272 
1273  RayTraversalData data;
1274
1275  // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1276
1277  while (!tstack.empty()) {
1278    data = tstack.top();
1279    tstack.pop();
1280
1281    if (!data.node->IsLeaf()) {
1282      // split the set of rays in two groups intersecting the
1283      // two subtrees
1284
1285      TraverseInternalNode(data, tstack);
1286     
1287    } else {
1288      // remove the ray from the leaf
1289      // find the ray in the leaf and swap it with the last ray...
1290      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
1291     
1292      if (!leaf->Mailed()) {
1293                leaf->Mail();
1294                if (affectedLeaves)
1295                  affectedLeaves->push_back(leaf);
1296       
1297                if (removeAllScheduledRays) {
1298                  int tail = leaf->rays.size()-1;
1299
1300                  for (int i=0; i < (int)(leaf->rays.size()); i++) {
1301                        if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1302                          // find a ray to replace it with
1303                          while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1304                                stat.removedRayRefs++;
1305                                leaf->rays[tail].mRay->Unref();
1306                                leaf->rays.pop_back();
1307                                tail--;
1308                          }
1309
1310                          if (tail < i)
1311                                break;
1312             
1313                          stat.removedRayRefs++;
1314                          leaf->rays[i].mRay->Unref();
1315                          leaf->rays[i] = leaf->rays[tail];
1316                          leaf->rays.pop_back();
1317                          tail--;
1318                        }
1319                  }
1320                }
1321      }
1322     
1323      if (!removeAllScheduledRays)
1324                for (int i=0; i < (int)leaf->rays.size(); i++) {
1325                  if (leaf->rays[i].mRay == ray) {
1326                        stat.removedRayRefs++;
1327                        ray->Unref();
1328                        leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1329                        leaf->rays.pop_back();
1330                        // check this ray again
1331                        break;
1332                  }
1333                }
1334     
1335    }
1336  }
1337
1338  if (ray->RefCount() != 0) {
1339    cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1340    exit(1);
1341  }
1342 
1343}
1344
1345
1346void
1347RssTree::AddRay(RssTreeNode::RayInfo &info)
1348{
1349
1350  stack<RayTraversalData> tstack;
1351 
1352  tstack.push(RayTraversalData(root, info));
1353 
1354  RayTraversalData data;
1355 
1356  while (!tstack.empty()) {
1357    data = tstack.top();
1358    tstack.pop();
1359
1360    if (!data.node->IsLeaf()) {
1361      TraverseInternalNode(data, tstack);
1362    } else {
1363      // remove the ray from the leaf
1364      // find the ray in the leaf and swap it with the last ray...
1365      RssTreeLeaf *leaf = (RssTreeLeaf *)data.node;
1366      leaf->AddRay(data.rayData);
1367      stat.addedRayRefs++;
1368    }
1369  }
1370}
1371
1372void
1373RssTree::TraverseInternalNode(
1374                                                          RayTraversalData &data,
1375                                                          stack<RayTraversalData> &tstack)
1376{
1377  RssTreeInterior *in =  (RssTreeInterior *) data.node;
1378 
1379  if (in->axis <= RssTreeNode::SPLIT_Z) {
1380
1381    // determine the side of this ray with respect to the plane
1382    int side = in->ComputeRaySide(data.rayData
1383                                                                  );
1384   
1385 
1386    if (side == 0) {
1387      if (data.rayData.mRay->HasPosDir(in->axis)) {
1388                tstack.push(RayTraversalData(in->back,
1389                                                                         data.rayData)
1390                                        );
1391                               
1392                tstack.push(RayTraversalData(in->front,
1393                                                                         data.rayData)
1394                                        );
1395       
1396      } else {
1397                tstack.push(RayTraversalData(in->back,
1398                                                                         data.rayData
1399                                                                         )
1400                                        );
1401                               
1402                tstack.push(RayTraversalData(in->front,
1403                                                                         data.rayData)
1404                                        );
1405                               
1406                               
1407      }
1408    } else
1409      if (side == 1)
1410                tstack.push(RayTraversalData(in->front, data.rayData));
1411      else
1412                tstack.push(RayTraversalData(in->back, data.rayData));
1413  }
1414  else {
1415    // directional split
1416        if (data.rayData.mRay->GetDirParametrization(in->axis - 3) <= in->position)
1417          tstack.push(RayTraversalData(in->back, data.rayData));
1418        else
1419          tstack.push(RayTraversalData(in->front, data.rayData));
1420  }
1421}
1422
1423
1424int
1425RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
1426{
1427  // first count all rays in the subtree
1428  // use mail 1 for this purpose
1429  stack<RssTreeNode *> tstack;
1430  int rayCount = 0;
1431  int totalRayCount = 0;
1432  int collapsedNodes = 0;
1433
1434#if DEBUG_COLLAPSE
1435  cout<<"Collapsing subtree"<<endl;
1436  cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1437  cout<<"depth="<<(int)sroot->depth<<endl;
1438#endif
1439
1440  //    tstat.collapsedSubtrees++;
1441  //    tstat.collapseDepths += (int)sroot->depth;
1442  //    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1443 
1444  tstack.push(sroot);
1445  VssRay::NewMail();
1446       
1447  while (!tstack.empty()) {
1448    collapsedNodes++;
1449    RssTreeNode *node = tstack.top();
1450    tstack.pop();
1451
1452    if (node->IsLeaf()) {
1453      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
1454      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1455                  ri != leaf->rays.end();
1456                  ri++) {
1457                               
1458                totalRayCount++;
1459                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
1460                  (*ri).mRay->Mail();
1461                  rayCount++;
1462                }
1463      }
1464    } else {
1465      tstack.push(((RssTreeInterior *)node)->back);
1466      tstack.push(((RssTreeInterior *)node)->front);
1467    }
1468  }
1469 
1470  VssRay::NewMail();
1471
1472  // create a new node that will hold the rays
1473  RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
1474  if (  newLeaf->parent )
1475    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1476 
1477
1478  tstack.push( sroot );
1479
1480  while (!tstack.empty()) {
1481
1482    RssTreeNode *node = tstack.top();
1483    tstack.pop();
1484
1485    if (node->IsLeaf()) {
1486      RssTreeLeaf *leaf = (RssTreeLeaf *) node;
1487     
1488      for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1489                  ri != leaf->rays.end();
1490                  ri++) {
1491                               
1492                // unref this ray from the old node
1493                               
1494                if ((*ri).mRay->IsActive()) {
1495                  (*ri).mRay->Unref();
1496                  if (!(*ri).mRay->Mailed()) {
1497                        (*ri).mRay->Mail();
1498                        newLeaf->AddRay(*ri);
1499                  }
1500                } else
1501                  (*ri).mRay->Unref();
1502                               
1503      }
1504    } else {
1505      tstack.push(((RssTreeInterior *)node)->back);
1506      tstack.push(((RssTreeInterior *)node)->front);
1507    }
1508  }
1509 
1510  // delete the node and all its children
1511  delete sroot;
1512 
1513  //   for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1514  //       ri != newLeaf->rays.end();
1515  //       ri++)
1516  //     (*ri).ray->UnMail(2);
1517
1518
1519#if DEBUG_COLLAPSE
1520  cout<<"Total memory before="<<GetMemUsage()<<endl;
1521#endif
1522
1523  stat.nodes -= collapsedNodes - 1;
1524  stat.rayRefs -= totalRayCount - rayCount;
1525 
1526#if DEBUG_COLLAPSE
1527  cout<<"collapsed nodes"<<collapsedNodes<<endl;
1528  cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1529  cout<<"Total memory after="<<GetMemUsage()<<endl;
1530  cout<<"================================"<<endl;
1531#endif
1532
1533  //  tstat.collapsedNodes += collapsedNodes;
1534  //  tstat.collapsedRays += totalRayCount - rayCount;
1535   
1536  return totalRayCount - rayCount;
1537}
1538
1539
1540int
1541RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1542{
1543  stack<RssTreeNode *> tstack;
1544  tstack.push(root);
1545
1546  Intersectable::NewMail();
1547  int pvsSize = 0;
1548       
1549  while (!tstack.empty()) {
1550    RssTreeNode *node = tstack.top();
1551    tstack.pop();
1552   
1553 
1554    if (node->IsLeaf()) {
1555          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1556          for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1557                  ri != leaf->rays.end();
1558                  ri++)
1559                if ((*ri).mRay->IsActive()) {
1560                  Intersectable *object;
1561#if BIDIRECTIONAL_RAY
1562                  object = (*ri).mRay->mOriginObject;
1563                  if (object && !object->Mailed()) {
1564                        pvsSize++;
1565                        object->Mail();
1566                  }
1567#endif
1568                  object = (*ri).mRay->mTerminationObject;
1569                  if (object && !object->Mailed()) {
1570                        pvsSize++;
1571                        object->Mail();
1572                  }
1573                }
1574        } else {
1575          RssTreeInterior *in = (RssTreeInterior *)node;
1576          if (in->axis < 3) {
1577                if (box.Max(in->axis) >= in->position )
1578                  tstack.push(in->front);
1579                               
1580                if (box.Min(in->axis) <= in->position )
1581                  tstack.push(in->back);
1582          } else {
1583                // both nodes for directional splits
1584                tstack.push(in->front);
1585                tstack.push(in->back);
1586          }
1587        }
1588  }
1589  return pvsSize;
1590}
1591
1592int
1593RssTree::CollectPvs(const AxisAlignedBox3 &box,
1594                                        ObjectContainer &pvs
1595                                        ) const
1596{
1597  stack<RssTreeNode *> tstack;
1598  tstack.push(root);
1599 
1600  Intersectable::NewMail();
1601 
1602  while (!tstack.empty()) {
1603    RssTreeNode *node = tstack.top();
1604    tstack.pop();
1605   
1606 
1607    if (node->IsLeaf()) {
1608          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1609          for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1610                  ri != leaf->rays.end();
1611                  ri++)
1612                if ((*ri).mRay->IsActive()) {
1613                  Intersectable *object;
1614                  object = (*ri).mRay->mTerminationObject;
1615                  if (object && !object->Mailed()) {
1616                        pvs.push_back(object);
1617                        object->Mail();
1618                  }
1619                }
1620        } else {
1621          RssTreeInterior *in = (RssTreeInterior *)node;
1622          if (in->axis < 3) {
1623                if (box.Max(in->axis) >= in->position )
1624                  tstack.push(in->front);
1625               
1626                if (box.Min(in->axis) <= in->position )
1627                  tstack.push(in->back);
1628          } else {
1629                // both nodes for directional splits
1630                tstack.push(in->front);
1631                tstack.push(in->back);
1632          }
1633        }
1634  }
1635  return pvs.size();
1636}
1637
1638void
1639RssTree::GetTreeStatistics(
1640                                                   float &avgPvsSize,
1641                                                   float &avgRays,
1642                                                   float &avgRayContribution,
1643                                                   float &avgPvsEntropy,
1644                                                   float &avgRayLengthEntropy,
1645                                                   float &avgImportance
1646                                                   )
1647{
1648  stack<RssTreeNode *> tstack;
1649  tstack.push(root);
1650       
1651  float sumPvsSize = 0.0f;
1652  float sumRayContribution = 0.0f;
1653  float sumImportance = 0.0f;
1654  float sumPvsEntropy = 0.0f;
1655  float sumRayLengthEntropy = 0.0f;
1656  float sumRays = 0.0f;
1657 
1658  int leaves = 0;
1659
1660  while (!tstack.empty()) {
1661    RssTreeNode *node = tstack.top();
1662    tstack.pop();
1663               
1664    if (node->IsLeaf()) {
1665          leaves++;
1666          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1667          leaf->UpdatePvsSize();
1668         
1669          sumPvsSize += leaf->GetPvsSize();
1670          sumRayContribution += leaf->GetAvgRayContribution();
1671          sumPvsEntropy += leaf->mPvsEntropy;
1672          sumRayLengthEntropy += leaf->mRayLengthEntropy;
1673          sumRays += leaf->rays.size();
1674         
1675          float imp = leaf->GetImportance();
1676         
1677          if (imp > 1.0f)
1678                cout<<"warning imp > 1.0f:"<<imp<<endl;
1679         
1680          sumImportance += imp;
1681         
1682
1683        } else {
1684          RssTreeInterior *in = (RssTreeInterior *)node;
1685          // both nodes for directional splits
1686          tstack.push(in->front);
1687          tstack.push(in->back);
1688        }
1689  }
1690 
1691  avgPvsSize = sumPvsSize/(float)leaves;
1692  avgRays = sumRays/(float)leaves;
1693  avgRayContribution = sumRayContribution/(float)leaves;
1694  avgPvsEntropy = sumPvsEntropy/(float)leaves;
1695  avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
1696  avgImportance = sumImportance/(float)leaves;
1697}
1698
1699
1700int
1701RssTree::GenerateRays(const float ratioPerLeaf,
1702                                          SimpleRayContainer &rays)
1703{
1704  stack<RssTreeNode *> tstack;
1705  tstack.push(root);
1706       
1707  while (!tstack.empty()) {
1708    RssTreeNode *node = tstack.top();
1709    tstack.pop();
1710               
1711    if (node->IsLeaf()) {
1712          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1713          float c = leaf->GetImportance();
1714          int num = (c*ratioPerLeaf + 0.5);
1715          //                    cout<<num<<" ";
1716
1717          for (int i=0; i < num; i++) {
1718                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1719                Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1720                Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1721                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1722                rays.push_back(SimpleRay(origin, direction));
1723          }
1724                       
1725        } else {
1726          RssTreeInterior *in = (RssTreeInterior *)node;
1727          // both nodes for directional splits
1728          tstack.push(in->front);
1729          tstack.push(in->back);
1730        }
1731  }
1732
1733  return rays.size();
1734}
1735
1736void
1737RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
1738{
1739  stack<RssTreeNode *> tstack;
1740  tstack.push(root);
1741       
1742  while (!tstack.empty()) {
1743    RssTreeNode *node = tstack.top();
1744    tstack.pop();
1745               
1746    if (node->IsLeaf()) {
1747          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1748          leaves.push_back(leaf);
1749        } else {
1750          RssTreeInterior *in = (RssTreeInterior *)node;
1751          // both nodes for directional splits
1752          tstack.push(in->front);
1753          tstack.push(in->back);
1754        }
1755  }
1756}
1757
1758bool
1759RssTree::ValidLeaf(RssTreeLeaf *leaf) const
1760{
1761  return true;
1762  //return leaf->rays.size() > termMinRays/4;
1763}
1764
1765
1766void
1767RssTree::PickEdgeRays(RssTreeLeaf *leaf,
1768                                          int &r1,
1769                                          int &r2
1770                                          )
1771{
1772  int nrays = leaf->rays.size();
1773
1774  if (nrays == 2) {
1775        r1 = 0;
1776        r2 = 1;
1777        return;
1778  }
1779 
1780  int tries = min(20, nrays);
1781
1782  while (--tries) {
1783        r1 = Random(nrays);
1784        r2 = Random(nrays);
1785        if (leaf->rays[r1].mRay->mTerminationObject != leaf->rays[r2].mRay->mTerminationObject)
1786          break;
1787  }
1788
1789  if (r1 == r2)
1790        r2 = (r1+1)%leaf->rays.size();
1791}
1792
1793
1794void
1795RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
1796                                                  const int numberOfRays,
1797                                                  SimpleRayContainer &rays)
1798{
1799
1800  int nrays = leaf->rays.size();
1801  for (int i=0; i < numberOfRays; i++) {
1802        bool useExtendedConvexCombination = (nrays >= 2) && (i > numberOfRays/2);
1803
1804               
1805        Vector3 origin, direction;
1806        // generate half of convex combination and half of random rays
1807        if (useExtendedConvexCombination) {
1808          // pickup 3 random rays
1809          int r1, r2;
1810          PickEdgeRays(leaf, r1, r2);
1811         
1812          Vector3 o1 = leaf->rays[r1].GetOrigin();
1813          Vector3 o2 = leaf->rays[r2].GetOrigin();
1814         
1815          const float overlap = 0.1f;
1816          float w1, w2;
1817          GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
1818          origin = w1*o1 + w2*o2;
1819          direction =
1820                w1*leaf->rays[r1].mRay->GetDir() +
1821                w2*leaf->rays[r2].mRay->GetDir();
1822          // shift the origin a little bit
1823          origin += direction*0.1f;
1824        } else {
1825          origin = GetBBox(leaf).GetRandomPoint();
1826          Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1827          direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1828        }
1829        //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1830        rays.push_back(SimpleRay(origin, direction));
1831  }
1832}
1833
1834int
1835RssTree::GenerateRays(const int numberOfRays,
1836                                          const int numberOfLeaves,
1837                                          SimpleRayContainer &rays)
1838{
1839
1840  vector<RssTreeLeaf *> leaves;
1841       
1842  CollectLeaves(leaves);
1843       
1844  sort(leaves.begin(),
1845           leaves.end(),
1846           GreaterContribution);
1847                         
1848
1849  float sumContrib = 0.0;
1850  int i;
1851  int k = 0;
1852  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
1853        if (ValidLeaf(leaves[i])) {
1854          float c = leaves[i]->GetImportance();
1855          sumContrib += c;
1856          //            cout<<"ray contrib "<<i<<" : "<<c<<endl;
1857          k++;
1858        }
1859                               
1860  float avgContrib = sumContrib/numberOfLeaves;
1861  float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
1862  k = 0;
1863  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
1864        if (ValidLeaf(leaves[i])) {
1865          k++;
1866          RssTreeLeaf *leaf = leaves[i];
1867          float c = leaf->GetImportance();
1868          int num = (c*ratioPerLeaf + 0.5);
1869          GenerateLeafRays(leaf, num, rays);
1870        }
1871
1872  return rays.size();
1873}
1874
1875
1876float
1877RssTree::GetAvgPvsSize()
1878{
1879  stack<RssTreeNode *> tstack;
1880  tstack.push(root);
1881
1882  int sumPvs = 0;
1883  int leaves = 0;
1884  while (!tstack.empty()) {
1885    RssTreeNode *node = tstack.top();
1886    tstack.pop();
1887               
1888    if (node->IsLeaf()) {
1889          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
1890          // update pvs size
1891          leaf->UpdatePvsSize();
1892          sumPvs += leaf->GetPvsSize();
1893          leaves++;
1894        } else {
1895          RssTreeInterior *in = (RssTreeInterior *)node;
1896          // both nodes for directional splits
1897          tstack.push(in->front);
1898          tstack.push(in->back);
1899        }
1900  }
1901
1902
1903  return sumPvs/(float)leaves;
1904}
1905
1906       
1907float
1908RssTreeLeaf::GetImportance()  const
1909{
1910
1911  if (1) {
1912        return GetAvgRayContribution();
1913        //      return GetPvsSize();
1914  } else {
1915        // return GetAvgRayContribution()*mEntropyImportance;
1916        //return GetAvgRayContribution();
1917        return mEntropyImportance;
1918  }
1919}
1920
1921
1922float
1923RssTreeLeaf::ComputePvsEntropy()
1924{
1925  int samples = 0;
1926  Intersectable::NewMail();
1927  // set all object as belonging to the fron pvs
1928  for(RssTreeNode::RayInfoContainer::iterator ri = rays.begin();
1929          ri != rays.end();
1930          ri++)
1931        if ((*ri).mRay->IsActive()) {
1932          Intersectable *object = (*ri).mRay->mTerminationObject;
1933          if (object) {
1934                if (!object->Mailed()) {
1935                  object->Mail();
1936                  object->mCounter = 1;
1937                } else
1938                  object->mCounter++;
1939                samples++;
1940          }
1941        }
1942 
1943  float entropy = 0.0f;
1944 
1945  if (samples > 1) {
1946        Intersectable::NewMail();
1947        for(RayInfoContainer::const_iterator ri = rays.begin();
1948                ri != rays.end();
1949                ri++)
1950          if ((*ri).mRay->IsActive()) {
1951                Intersectable *object = (*ri).mRay->mTerminationObject;
1952                if (object) {
1953                  if (!object->Mailed()) {
1954                        object->Mail();
1955                        float p = object->mCounter/(float)samples;
1956                        entropy -= p*log(p);
1957                  }
1958                }
1959          }
1960        entropy = entropy/log((float)samples);
1961  }
1962  else
1963        entropy = 1.0f;
1964 
1965  return entropy;
1966}
1967
1968float
1969RssTreeLeaf::ComputeRayLengthEntropy()
1970{
1971  // get sum of all ray lengths
1972  // consider only passing rays or originating rays
1973  float sum = 0.0f;
1974  int samples = 0;
1975  int i=0;
1976  for(RayInfoContainer::const_iterator ri = rays.begin();
1977      ri != rays.end();
1978      ri++)
1979        if ((*ri).mRay->IsActive()) {
1980//        float s;
1981//        if (i == 0)
1982//              s = 200;
1983//        else
1984//              s = 100;
1985//        i++;
1986         
1987          sum += (*ri).mRay->GetSize();
1988          samples++;
1989        }
1990 
1991 
1992  float entropy = 0.0f;
1993 
1994  if (samples > 1) {
1995        i = 0;
1996        for(RayInfoContainer::const_iterator ri = rays.begin();
1997                ri != rays.end();
1998                ri++)
1999          if ((*ri).mRay->IsActive()) {
2000//              float s;
2001//              if (i==0)
2002//                s = 200;
2003//              else
2004//                s = 100;
2005//              i++;
2006//              float p = s/sum;
2007                float p = (*ri).mRay->GetSize()/sum;
2008                entropy -= p*log(p);
2009          }
2010        entropy = entropy/log((float)samples);
2011  } else
2012        entropy = 1.0f;
2013 
2014  return entropy;
2015}
2016 
2017
2018
2019void
2020RssTreeLeaf::ComputeEntropyImportance()
2021{
2022  mPvsEntropy = ComputePvsEntropy();
2023  mRayLengthEntropy = ComputeRayLengthEntropy();
2024 
2025  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
2026  mEntropyImportance = 1.0f - ComputePvsEntropy();
2027 
2028  //  cout<<"ei="<<mEntropyImportance<<" ";
2029}
2030
2031
2032AxisAlignedBox3
2033RssTree::GetShrankedBBox(const RssTreeNode *node) const
2034{
2035  if (node->parent == NULL)
2036        return bbox;
2037 
2038  if (!node->IsLeaf())
2039        return ((RssTreeInterior *)node)->bbox;
2040
2041  // evaluate bounding box from the ray origins
2042  AxisAlignedBox3 box;
2043  box.Initialize();
2044  RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2045  for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2046          ri != leaf->rays.end();
2047          ri++)
2048        if ((*ri).mRay->IsActive()) {
2049          box.Include((*ri).GetOrigin());
2050        }
2051  return box;
2052}
2053
2054
2055int
2056RssTree::CollectRays(VssRayContainer &rays,
2057                                         const int number)
2058{
2059  VssRayContainer allRays;
2060  CollectRays(allRays);
2061 
2062  int desired = min(number, (int)allRays.size());
2063  float prob = desired/(float)allRays.size();
2064  while (rays.size() < desired) {
2065        VssRayContainer::const_iterator it = allRays.begin();
2066        for (; it != allRays.end() && rays.size() < desired; it++) {
2067          if (Random(1.0f) < prob)
2068                rays.push_back(*it);
2069        }
2070  }
2071  return rays.size();
2072}
2073
2074
2075int
2076RssTree::CollectRays(VssRayContainer &rays
2077                                         )
2078{
2079  VssRay::NewMail();
2080
2081  stack<RssTreeNode *> tstack;
2082  tstack.push(root);
2083 
2084  while (!tstack.empty()) {
2085    RssTreeNode *node = tstack.top();
2086    tstack.pop();
2087       
2088    if (node->IsLeaf()) {
2089          RssTreeLeaf *leaf = (RssTreeLeaf *)node;
2090          // update pvs size
2091          RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
2092          for (;it != leaf->rays.end(); ++it)
2093                if (!(*it).mRay->Mailed()) {
2094                  (*it).mRay->Mail();
2095                  rays.push_back((*it).mRay);
2096                }
2097        } else {
2098          RssTreeInterior *in = (RssTreeInterior *)node;
2099          // both nodes for directional splits
2100          tstack.push(in->front);
2101          tstack.push(in->back);
2102        }
2103  }
2104 
2105  return rays.size();
2106}
Note: See TracBrowser for help on using the repository browser.