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

Revision 463, 46.8 KB checked in by bittner, 19 years ago (diff)

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