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

Revision 1883, 47.3 KB checked in by bittner, 18 years ago (diff)

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