source: trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp @ 473

Revision 473, 46.7 KB checked in by mattausch, 19 years ago (diff)

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