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

Revision 447, 46.5 KB checked in by bittner, 19 years ago (diff)

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