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

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

test version
build script for testing. code is set uo for testing also
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 minPvs )\n"<<
180    minPvsNodes*100/(double)Leaves()<<endl;
181
182  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minRays )\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 = (float)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.