source: GTP/trunk/Lib/Vis/Preprocessing/src/VssTree.cpp @ 1824

Revision 1824, 47.2 KB checked in by bittner, 18 years ago (diff)

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